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1 . GENERAL 

The work completed under the Delivery order contract DO-41 comprised the final part of a 

three year study of the sensitivity of low gravity experiments to residual acceleration experienced 
aboard a spacecraft in low earth orbit. The sensitivity analyses were carried out using “generic 
models” which were representative of the types of experiment that would take place on Space 
Station. This report summarizes all the work carried out over the period 7/5/89-10/4/90 . 

The following papers have been published, accepted for publication or are in preparation 
for submission for publication: 

1. A. Nadarajah, F. Rosenberger and J.I.D. Alexander, Modelling the solution growth of 
triglycine sulfate in low gravity, J. Crystal Growth 104 (1990) 218-232. 

2. J. I. D. Alexander and F. Rosenberger, Bridgman crystal growth in low gravity: A scaling 
analysis, to appear in Low Gravity Fluid Dynamics and Transport Phenomena, J. Koster 
and R. Sani (eds.) AIAA, (1990). 

3. J.I. D. Alexander and J. Ouazzani, A Pseudospectral Collocation Method Applied to the 
Problem of Convective-Diffusive Transport in Fluids Subject to Residual Acceleration, 
Proc. 6th International Conference on Numerical Methods in Laminar and Turbulent Flow, 
C. Taylor, P. Gresho, R. L. Sani and J. Hauser (eds.) Pineridge, Swansea (1989) 1035. 

4. J. I. D. Alexander, Low-Gravity Experiment Sensitivity to Residual Acceleration: A 
Review, Microgravity Science and Technology 3, (1990) 52-68. 

5. F. Rosenberger, J.I.D. Alexander, A. Nadarajah and J. Ouazzani, Influence of Residual 
Gravity on Crystal Growth Processes, to appear Microgravity Science and Technology 
(1990). 

6. J.I.D. Alexander, S. Amiroudine, J. Ouazzani and F.Rosenberger F.: Analysis of the Low 
Gravity Tolerance of the Bridgman-Stockbarger Technique. II: Transient and Periodic 
Accelerations, to be submitted to the J. Crystal Growth (1990). 

7. A. Nadarajah, F. Rosenberger and J.I.D. Alexander, Modelling Crystal Growth by 
Physical Vapor Transport in Closed Ampoules, to be submitted to the J. Crystal Growth 
(1990). 

In addition to the above publications, the results of our work have been presented at the 
following conferences 

1 . “Commercial Numerical Codes: To Use or Not to Use, Is This The Question?” presented 
by J.I.D Alexander at the Microgravity Fluids Workshop, Westlake Holiday Inn, Cleveland Ohio, 
August 7-9, 1990. 

2. “Sensitivity of Crystal Growth Experiments to Residual Accelerations, presented by J.I.D. 
Alexander, the Gordon Conference on Gravitational Effects in Materials and Processes”, July 30 - 
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August 4, 1989, Plymouth State College, Plymouth, New Hampshire. 

3. “Low Gravity Accelerations: Influence on Dopant Redistribution During Directional 
Solidification”, presented by J.I.D. Alexander at the Third International Symposium on 
Experimental Methods for Microgravity Materials Science Research, COSPAR XXVIII Plenary 
Meeting, June 25-July 6, 1990. 

4. “Modelling the Solution Growth of TGS Crystals in Low Gravity”, presented by J.I.D 
Alexander, at the Third International Symposium on Experimental Methods for Microgravity 
Materials Science Research, COSPAR XXVIII Plenary Meeting, June 25-July 6, 1990. 

5 . “Modelling the Solution Growth of TGS Crystals in Low Gravity”, presented by A. 
Nadarajah at the Eighth American Conference on Crystal Growth, Vail, Colorado July (1990). 

6. “Influence of Residual Gravity on Crystal Growth Processes,” presented By F. 
Rosenberger at the First International Microgravity Congress, Bremen September (1990). 

The results of the work carried out under DO-41 are presented in the following sections. 
The work has mostly been concerned with the effects of time-dependent accelerations during 
directional solidification by the Bridgman technique, steady acceleration during crystal growth by 
physical vapor transport (PVT) and the completion of a study of steady and time-dependent 
acceleration effects during TGS solution growth. In section 2 we review work related to the 
analysis of the sensitivity of residual acceleration. In section 3 we present the latest results for the 
analysis of the sensitivity of the Bridgman technique to time-dependent residual acceleration. The 
results include the effect of oscillatory accelerations as well as acceleration time-series synthesized 
from measured disturbances on spacelab-3. In section 4 we present the results of our analysis of 
the process of crystal growth by Physical Vapor Transport. The results of the study of the TGS 
solution growth system are given in the appendix (section 5.2) 


2. EXPERIMENT SENSITIVITY TO RESIDUAL ACCELERATION: A REVIEW 
2.1. Introduction 

It has been recognized for some time that the low-gravity acceleration environment associated 
with a spacecraft in low earth orbit offers an opportunity to study certain physical processes which 
are difficult to investigate under the gravitational acceleration experienced at the earth's surface. 
Our experience with such opportunities has led us to realize that the residual acceleration 
environment on board an orbiting spacecraft is not as low or as steady as would be desired for 
certain classes of experiments . The sources of the residual acceleration include [1-5] crew 
motions, mechanical vibrations (pumps, motors, excitations of natural frequencies of spacecraft 
structures), spacecraft maneuvers and basic attitude motion, atmospheric drag and the earth's 
gravity gradient. The accelerations are characterized by temporal variations in both magnitude and 
orientation. Such disturbances will, in particular, affect those experiments susceptible to buoyancy 
induced fluid motion [6]. Indeed, recent studies [7,9] indicate that transient disturbances can have 
undesirable long term effects. The analysis of the sensitivity of any proposed flight experiment is 
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necessary for a variety of reasons. These include the need to optimize the use of the limited time 
available for flight experiments, the interpretation of experimental results and the determination of 
tolerable acceleration levels to be used in planning for NASA's Space Station [3]. 

Examples of experiments that are conducted in the low-gravity environment include critical 
point studies [10], crystal growth [11-18], and diffusion experiments [19,20]. Since these 
experiments involve density gradients in the fluid phase, they are inherently sensitive to the effects 
of buoyancy-driven fluid motion [21-25]. For certain types of systems it is known [7,8,26-37] that 
the associated physical processes (mass transfer, heat transfer, convection, oscillation and 
distortion of free surfaces) are sensitive to time-dependent accelerations. For example, the transfer 
of mass and heat in fluid systems can be significantly affected by oscillatory flows [28-37]. The 
effect of oscillation enhanced transport conditions on the local variation in composition during 
crystal growth is not well characterized for most of the systems relevant to crystal growth under 
microgravity conditions. Order of magnitude analyses [7,8,38] suggest that for certain 
combinations of physical properties and growth conditions, oscillations in the residual acceleration 
may adversely affect the mass transport conditions. For instance, the DMOS experiment [39] on 
STS 51-D showed extensive evidence of convective mixing of liquids. The degree of mixing is 
greater than can be attributed to the quasi-steady low-gravity component, but can be accounted for 
when oscillations in effective gravity are included in the description of the transport conditions. 

It is possible that in future space laboratories some experiments will require isolation from 
vibration. The design of efficient isolation systems, however, requires determination of the 
tolerable acceleration levels (amplitudes and frequencies) for given experiments. The results of 
sensitivity analyses can be used in the development of vibration isolation systems to identify those 
experiments that will need isolation and by supplying data concerning frequency dependence, the 
effect of transients etc.. The system can then be designed to filter out those bandwidths to which 
the experiment is predicted to respond adversely. 

Low-gravity experiments which involve liquids with free surfaces are most likely to be 
sensitive to vibrations. Experiments of this type include studies of liquid bridges, equilibrium 
shapes of drops and bubbles [40-43] and thermocapillary flow experiments [43-46]. Evidence for 
the oscillation of liquid zone shapes was found on the D-l mission. Long liquid bridges were more 
sensitive to residual gravity, and exhibited random oscillations in zone shape [40]. 

It should be recognized that the sensitivity of a given process to the overall low-gravity 
environment will depend on one or more of the following factors: 

a) magnitude and direction of accelerations 

b) system geometry 

c) boundary conditions (eg. insulated or conducting walls, solute sources and sinks), and 

d) physical properties of the participating materials (viscosity, thermal and solute 
diffusivities). 

Any analysis of the sensitivity to oscillatory accelerations and transient disturbances characteristic 


3 



of the anticipated environment aboard Space Station, or any other spacecraft, must take these 
factors into account. In general, a full mathematical characterization may not be practical. 
However, identification of sensitive experiments in combination with a judicious choice of 
simplified numerical models can be used to choose operating parameters that will optimize the use 
of the low-gravity environment. With this in mind, we review work related to the effects of steady 
and time-dependent residual accelerations, and associated work concerning transport in oscillating 
flows. 

In section 2 we describe the components of the low-gravity environment. In section 3 
previous work involving order of magnitude (OM) estimates is discussed. The application of OM 
analyses to acceleration-sensitivity determination is discussed in 3.1. The results of simple 
analytical models are summarized in section 3.2. Work related to linear stability in the presence of 
modulated gravity is briefly discussed in section 3.3. The results of direct numerical simulations 
are reviewed in section 3.4. In section 4., the results and utility of residual acceleration analyses 
are summarized and discussed. 

2.2. The microgravity environment 

The low-gravity environment of a spacecraft has been referred to as a "zero gravity" 
environment. This is because any object within the craft is subject to roughly the same acceleration 
due to the earth's gravitational field. Thus, a free object may move along approximately the same 
orbital path as the spacecraft. As a consequence, to an observer in the space craft objects appear to 
behave as if no gravity were present. However, any body capable of motion relative to the 
spacecraft will experience an acceleration relative to the mass center of the spacecraft. This relative 
acceleration arises from several sources which include the gravity gradient of the earth, 
atmospheric drag and attitude motions, as well as machineiy vibrations and crew activities. 

The relative acceleration that is associated with the earth's gravity gradient is defined as 
follows. The mass center of the spacecraft is subject to a force F(r 0 ) due to the gravitational 
attraction of the earth, where r 0 is the position of the mass center of the craft with respect to the 
mass center of the earth. A particle at a position r within the spacecraft is subject to a force F(r). 
If, as a first approximation, we ignore the interaction between the spacecraft and the mass of the 
particle, and if the distance |r - r 0 1 is small compared to the semi-major axis of the spacecraft orbit, 
the difference between the forces applied at r and r 0 define the gravity gradient G as follows [2-6] 

F(r) - F(r 0 ) = Gz, z = r - r 0 . (1) 

If the spacecraft maintains a fixed orientation, for example with respect to the sun, there is no 
rotation of the spacecraft frame relative to the geocentric frame. In this type of orbit both the 
magnitude and the orientation of the residual acceleration vector change as a function of time and 
there will be no steady residual acceleration. 

In addition to the gravity gradient acceleration the atmospheric drag force on the spacecraft 
can be significant [2-6], despite the fact that the atmospheric density at typical shuttle altitudes is 
only on the order of 10" 12 kg m~ 3 . Atmospheric drag causes a slow inward spiral of the spacecraft. 
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To an observer in the spacecraft the effect of atmospheric drag is to produce a relative acceleration 
that is equal in magnitude but opposite in sign to the atmospheric drag force. Table 1 gives 
atmospheric drag acceleration estimates and relative acceleration related to the gravity gradient and 
centrifugal force in |ig per meter from the mass center of a spaceraft that is in a circular orbit with a 
gravity gradient stabilized attitude. Xi is parallel to the local vertical and X 2 is perpendicular to the 
orbital plane. For this attitude the acceleration points away from the mass center along Xi and 
toward the mass center along X 2 . 

If the spacecraft does not maintain such a fixed orientation then, in addition to (1), 
centrifugal, Coriolis and Euler accelerations will become apparent in the spacecraft reference frame. 
They arise in conjuction with spacecraft attitude motions and depend on the nature of the rotation of 
the spacecraft frame of reference relative to a fixed geocentric frame. If a spacecraft in a quasi- 
circular orbit continuously rotates relative to the fixed geocentric frame such that a given direction 
in the spacecraft frame is always oriented parallel to the craft’s position vector r 0 , a centrifugal 
acceleration will augment the gravity gradient in the direction parallel to the position vector, and 
cancel the component tangent to the flight path. Under these conditions there will be a steady 
component to the residual acceleration vector. In practice this situation arises, for example, in the 
so-called gravity gradient stabilized attitude. 

Euler accelerations arise when the rate of rotation of the spacecraft frame changes as a 
function of time, as it would in a non-circular orbit, or as the semi-major axis of the orbit changes 
as the craft slowly spirals inward due to atmospheric drag. Coriolis accelerations will also be 
apparent but are proportional to the velocity of an object relative to the spacecraft. Thus, unless 
relative velocities are large, the Coriolis effect will be insignificant [6]. 

In addition to the above, there are a variety of residual acceleration components that occur 
over a broad range of amplitudes and frequencies. The nature and sources of these accelerations 
have been discussed in [1-3]. Figures 1 and 2 are examples of residual accelerations measured on 
orbit. 


2.3. Analyses of residual acceleration effects 

As discussed in the introduction, many experiments are undertaken in orbiting spacecraft in 
order to reduce or eliminate unwanted effects of buoyancy driven motions in fluids. Sources of 
buoyancy in fluids include discrete interfaces, as they occur, for example between immiscible 
fluids, solid particles and fluids, and bubbles and fluids with different densities, other buoyancy 
sources are internal density gradients that arise as a consequence of gradients in temperature and 
concentration. The nature of the response of a bulk fluid to the effective body force associated 
with a low-gravity environment will depend, in part, on the relative orientation of the residual 
acceleration vector and the density gradient. In particular, if the integral of the effective body force 
pb (where p is the density and b the effective acceleration) acting on any closed circuit within the 
fluid is non-zero, motion will immediately ensue. That is, if 
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( 2 ) 


j dc-pb = j dA-curl pb = j dA-[grad pAb + p curl b] * 0, 

motion will occur. Note also, that even if there is no density gradient motion will occur if 

dA-curl pb *0. (3) 


» 

s 


This includes the case of "precessional stirring" [48]. This could occur, for example, when a 
spacecraft orbiting in a radial attitude (i.e. it continuously rotates about an axis perpendicular to the 
orbital plane) underwent an additional rotation about some axis other than the basic orbital axis. 

There have been three appraoches to the determination of experiment sensitivity to the 
residual acceleration environment. 

a) order of magnitude analyses [38,49-52], 

b) analytical models [5,6,39,53-59], and 

c) direct numerical simulation [9,3,60-73]. 

In addition to work directly related to microgravity, there is extensive literature on mass and 
heat transport in oscillatory flows [26-36], the effect of gravity modulation on the classical Benard 
problem [74-77] and the linear stability of planar fluid surfaces [78,79]. 

2.3.1. Order of magnitude analyses 

An order of magnitude (OM) or dimensional analysis of a physical process involves the 
examination of the system of equations that are assumed to govern the process. For fluid physics 
and material science experiments, the relevant physical processes are governed by equations 
describing the transport of heat, mass and momentum in single and multi-component fluid 
systems. The usual approach to OM analyses [51] involves the definition of reference quantities 
which appropriately characterize the physical system (i.e. reference scales for velocity, time, 
length, forces, etc.). Application of these scales to the (dimensional) governing equations, 
followed by the definition of characteristic dimensionless groups admits a comparison of the order 
of magnitude of every term in each equation. There are several comprehensive discussions of OM 
analyses applied to low-gravity situations (eg. [38,49-51]). 

The advantage of the OM approach is that a great deal of information pertaining to the 
sensitivity of a given experiment can be obtained with little computational effort. The 
disadvantages are related to the fact that the choice of characteristic reference quantities for length, 
velocity, time etc. are not, in general, known a priori [80]. Thus, in cases for which the chosen 
characteristic scale is not (for a given set of conditions) representative of the actual process, the 
results of the analysis may be in error, often by several orders of magnitude. In addition, 
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application of OM analyses to the problem of g-sensitivity has invariably involved implicit 
linearization and often poorly represents the multi-dimenional nature of the physical process. This 
can also lead to incorrect OM predictions. Finally, analyses are generally restricted to the 
examination of a single component disturbance (e.g. ~ sin(cot) ). As a result, they may not properly 
indicate the response of the system to the typically complex multi-frequency disturbances 
characteristic of the spacecraft environment. It is nonetheless useful to examine the general trends 
predicted by such analyses since they can provide at least qualitative guidance for more detailed 
analytical and numerical studies. 

The determination of tolerable g-levels using OM analyses is based upon estimations of the 
response of a system to either steady accelerations or time-dependent disturbances of the form 

g(t) = a cos(cot), (4) 

where co is the angular frequency of the disturbance. The susceptibility of an experimental system 
to such disturbances is defined via the magnitude of a particular response (say a temperature, 
compositional or velocity fluctuation) which must not be exceeded in order to bring the experiment 
to a successful conclusion. The most obvious trend predicted by analyses to date is shown in 
Figs. 3-5. The curves depict the maximum tolerable residual acceleration magnitude as a function 
of frequency, /= co/27t, for given experiments and have the form 


gmax — F(CO,n crit , pj) , \ ) 

where n cr it is the magnitude of the maximum allowable response, and the pi, i=l,N represent the 
N material properties of the system (such as viscosity, thermal and species diffusivities). In 
general, below 10" 2 to 10'3 Hz the response is expected to be either independent or weakly 
dependent on frequency. At higher frequencies the analyses indicate that the allowable residual 
acceleration magnitude increases as approximately the square of the frequency. This behavior can 
be understood upon consideration of the following OM estimates [38] (used to obtain Fig. 3) for 
velocity, temperature and concentration changes associated with the acceleration (4): 

V ~ (6) 

(co 2 + v 2 /L 4 ) 1/2 


5T~ 


Vg 


AT/L 

(CD 2 + K 2 /L 4 ) 1/2 ’ 


(7) 


8C 


V ac/l 

V + D 2 /L 4 ) 1/2 ’ 


( 8 ) 
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where g* is the residual acceleration magnitude, p is the mass density, and Ap/p is the relative 
change in density due to temperature and/or composition; AT/L, AC/L and L are, respectively, a 
characteristic temperature gradient, concentration gradient and length; k is the thermal diffusivity. 
The sensitivity limits are obtained by using either 8T, 8C or V g as a sensitivity parameter and 
setting g* equal to the allowable residual acceleration gmax- 

For situations which involve mixed thermocapillary convection and buoyancy-driven 
convection the governing dimensionless groups are [51]; the Bond number, B 0 = pg*L 2 /y, the 
Grashof number Gr = AT |J3t I g*L 3 /v 2 > the Prandtl number Pr = v /k, and the Surface Reynolds 
number Rs = I y T lATL/(o.v. Here I y T I is the rate of change of surface tension y with temperature, 
p T is the coefficient of thermal expansion, v is the kinematic viscosity, ji is the shear viscosity 
These groups respectively represent the ratios of buoyancy force to surface force, buoyancy force 
to inertial force, thermal and momentum timescales, and surface force to inertial force. 

The relative importance of gravity and thermocapillary forces can then be estimated upon 
considering the relative magnitudes of the velocity due to buoyancy and the velocity due to 
thermocapillary effects [50,51]. The latter is given by 


y* _ NAT 


(9) 


For thermocapillary flow to dominate, the ratio of V g to V* must be less than one. For this to 
occur, the maximum magnitude of the acceleration must satisfy: 


VpT 


GO) 


Equation (10) can be re-interpreted in terms of the ratio of the Grashof number, Gr, and the 
surface Reynolds number Rs. The condition for buoyancy and thermocapillary forces to be the 
same magnitude is 

Gr = Rs(St 2 + 1) 1/2 , (11) 

where St= coL 2 /v is the Strouhal number. This condition is illustrated in Fig. 6. Three regimes of 
interest are identified. These regions are defined by the relative values of the Grashof, Strouhal and 
surface Reynolds numbers. In region I, the Strouhal number St is less than one, and the condition 
that thermocapillary forces dominate is equivalent to that for a steady flow. In region II, St is 
greater than one. In this region, for a fixed value of Rs, the value of Gr required to give buoyancy 
forces equal weight to thermocapillary forces increases with increasing values of the Strouhal 
number. This essentially reflects the fact that the characteristic time for the fluid response greatly 
exceeds the period of the disturbance. Thus, thermocapillary forces dominate at higher values than 
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in region I. In region in buoyancy forces predominate over thermocapillary forces. For situations 
here a boundary layer scaling is appropriate i.e Rs » 1 [50]. The relative importance of buoyancy 
can be estimated via the ratio 


Gr = Rs 5/3 (St 2 + 1) 1/2 . (12) 

Figure 7 depicts the relationship between the estimated tolerable g-level and the frequency of 
the acceleration for experiments with the physical properties given in Table 2. These figures 
represent estimates based on two extreme flow regimes, penetrative or "viscous" flow (Fig. 7a) 
and boundary layer flow (Fig. 7b) and may be respectively considered as "worst" and "best" case 
estimates for the given sensitivity parameters. The curves have been determined for a sensitivity 
criterion which requires the tolerable acceleration levels to be such that the buoyancy forces should 
not exceed 10% of the thermocapillary forces (if the requirement had been 1%, the tolerable 
accelerations would be an order of magnitude smaller). In terms of acceleration levels predicted for 
the space station or measured on past missions the practical sensitive ranges are restricted to 
disturbances with frequencies less than 10 Hz. The majority of the estimates suggest that ,for these 
experimental parameters, periodic vibrations will generally not lead to significant buoyancy effects. 

The above examples estimate the effects of simple single component accelerations. In reality, 
low-gravity disturbances tend to be associated with more than one frequency. Given that a system 
may respond to a multi-frequency disturbance in an "additive" way, tolerance curves such as these 
may underestimate the response to a given low-gravity disturbance. For example, models of the 
DMOS experiment show that mass transport has an additive response to the (multi-frequency) 
residual acceleration [39]. Post flight analysis of the experiment results has demonstrated that the 
amount of mixing observed between organic liquids can be explained by the additive response of 
the system to a multicomponent disturbance. 

In many cases an a priori choice of length or time scales may not be obvious. For example, 
in a study of dopant (solute) uniformity in directionally solidified crystals, the OM estimates of 
Rouzaud et al. [52] and Camel and Favier [83] are in agreement with the direct numerical 
simulations of Chang and Brown [84] for a Schmidt number (ratio of the melt's kinematic 
viscosity to dopant diffusivity) of fifty but overestimate the amount of radial segregation for a 
Schmidt number of ten. 

The results of an order of magnitude or scaling analysis have been compared in [80] with 
those of numerical simulations of the effects of steady residual acceleration on compositional 
nonuniformity in directionally solidified crystals. The basic model consisted of a Bridgman-type 
system with an ampoule translated between hot and cold zones separated by a thermal barrier. In 
order to apply the Camel-Favier technique [52,83] to the numerical simualtion, it was necessary to 
evaluate a proportionality constant, x 0 2 . This relates the maximum convective velocity obtained in 
the numerical simulations with the Grashof number. The value of this constant depended on the 
orientation of the acceleration, and on the chosen length scale. The adiabatic zone length L a yielded 
results which best showed the trends predicted by the Camel-Favier approach. Table 3 lists values 
of x 0 associated with different choices of reference length. Having calculated x 0 we then 
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determined the transport regime (see references [52,80,83] for details) by graphing, in Fig. 8, 
GrSc vs Pe g for the case of the residual acceleration parallel to the interface. Here Sc (Schmidt 
number) gives the ratio of kinematic viscosity to thermal diffusivity, and Pe g = VL/D, is a ratio of 
the translation rate to the characteristic diffusive speed where L is the characteristic diffusion length 
and D the solute diffusivity. The lateral composition non-uniformities in the crystal obtained from 
numerical simulations [9] are also given in Fig. 8. A comparison between computed and predicted 
nonuniformity values showed that at high values of GrSc (region Ic in Fig. 8), the order of 
magnitude approach meets with some success but is less faithful at low values of GrSc (for 
example the transition between regions lb and II). Fig. 8 also contains two cases for which Pe g is 
reduced at fixed GrSc. According to the estimates the lower GrSc case should have yielded a 
higher compositional nonuniformity as the growth rate was reduced, i.e. the estimates predict that 
the system moves closer to the convective diffusive transition. This was not reflected in the 
nonuniformity obtained from the numercial simulations. Evaluation of the numerical results 
showed that the increase in the characteristic diffusion length scale (associated with reduction in 
growth rate) led to an order of magnitude decrease in the solute gradient at the interface. Hence, 
because the system-wide variation in composition is small, the interfacial nonuniformity is 
correspondingly small. 

The above results indicate that scaling arguments alone can at best be expected to yield order 
of magnitude accuracy. The limitations of scaling arguments are hardly surprising if one considers 
the neglect of the multi-dimensionality (boundary conditions, flow structure, etc.) of the physical 
situation which is inherent in such scalar descriptions. Of course, prior knowledge of the system 
behavior can be used to locally improve the accuracy of the scaling, (eg. Kimura and Bejan [85]). 
Similarly, if a system locally exhibits one-dimensional behavior, such as boundary layer flow, then 
the appropriate scaling for the formal reduction of the transport equations can be used effectively to 
estimate the local system response. In general, however, our understanding of the response of 
systems to residual acceleration can be furthered only by an approach which uses scaling, 
mathematical modelling and, naturally, the results of experiments in a complementary fashion. 

2.3.2. Analytical models 

A few attempts have been made to assess the effects of time-dependent residual acceleration 
using analytical models [53-59]. The first attempt to model the effects of disturbances on heat 
transfer between a fluid and a solid bounding medium was carried out by Gebhardt [53]. Under 
the assumption that the fluid and its container would be subject to a "sequence of abrupt relative 
displacements spaced by a time interval x c , ...", Gebhardt took x c to be a random distributed 
variable. He showed that as x m (the most likely value of x c ) decreases, the resulting heat transfer 
was greater than predicted by conduction in the absence of a disturbance. In particular, it was 
found that for x m < d 2 /K (where d is a characteristic length and k is the thermal diffusivity) the heat 
transfer increased rapidly as Xm decreased. The main limitation of this model is that fluid 
convection (and thus fluid properties such as kinematic viscosity) does not enter the analysis. It is 
well known that the response of heat transfer to convection varies according to the Prandtl number. 
In particular, for low-gravity flows (as we shall see later) the temperature field is less sensitive to 
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convection for low Prandtl number fluids (Pr < 1). 

Approximate solutions for transient convection in a cylinder with an azimuthal temperature 
variation have been obtained by Dressier [54]. Residual accelerations representing the motion of an 
astronaut and a transient rotation of the spacecraft were imposed on the system. Linear 
accelerations were imposed perpendicular to the cylinder axis. For kinematic viscosities of 10" 2 
cm 2 s'l, cylinders of radius 1 and 2 cm with a maximum temperature difference of 95 K across the 
diameter were examined. The maximum fluid velocities consequent to the astronaut motion 
(modelled by a sequence of two step functions separated by 2 seconds with 4 x 10‘ 4 g magnitude, 
0.5 second duration and zero mean, oriented perpendicular to the temperature gradient) were 12 
pm s' 1 for both the large and small diameter cylinders (see Fig. 9). For the larger cylinder the 
second pulse reduced the velocity from 12 pm s" 1 to 0 by the end of the pulse. The effect of the 
second pulse was to change the velocity in the small container from 12 to -3 pm s' 1 by the end of 
the pulse. The velocity decayed asymptotically to zero reaching a value close to zero after 15 
seconds. The reason the response was different for the two containers can be explained simply in 
terms of the momentum decay times (ti /2 ~ R 2 /v) for each container. The velocity in the larger 
container had hardly begun to decay when the second impulse arrived. The shorter decay time 
associated with the smaller container meant that the impulse acted on a smaller magnitude velocity 
which explains the change in direction of the motion. While the disturbance had a mean value of 
zero the fluid particles attained a finite displacement (indicated by the value x of the angular 
displacement in Fig. 6). Similarly, a disturbance corresponding to a simple spacecraft maneuver 
was examined. The simulated maneuver corresponded to a 90° rotation of the spacecraft at a rate of 
1° s' 1 . The maximum velocities associated with the centrifugal acceleration were, respectively, 40 
and 35 pm s' 1 for the large and small diameter cylinders. 

Experiments involving free surfaces of liquids have a high probability of responding to 
residual accelerations. Indeed, unexpected results have been attributed to residual accelerations in 
several cases [40,86-88]. Particular examples of such experiments involve liquid (float) zones. 
Three types of liquid zones are typically the subject of low-gravity experiments [89]: isothermal 
and non-isothermal zones suspended between inert solids (i.e. liquid bridges), and non-isothermal 
zones suspended between a feed rod and a growing crystal. The dimensionless group associated 
with the shape of isothermal static liquid zones is the Bond number, B 0 , defined earlier. It 
expresses the relative importance of gravitational and surface forces. A comprehensive description 
of stable zone shapes in zero gravity, and as a function of B 0 , is given by Martinez et al. [41]. 
Fig. 10 displays their results for the stability of static zones with a steady acceleration parallel to the 
zone axis. In order to relate the Bond number to a given acceleration level they give the value of 
B 0 = 1.4 for a 1 liter volume of water subject to an acceleration of 1 cm s' 2 . For B 0 = 0, the zone 
volume V=V*= tcR 2 L and contact angle (j)= 90° the maximum stable length is 27tR (the Rayleigh 
limit). Coriell et al. [90] investigated the effect of B 0 on the maximum length of such zones and 
compared theoretical predictions with experimental results (see Fig. 1 1). The actual length of the 
equilibrium zone is determined by the volume, the contact angle and the static Bond number. 
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Estimates of the effect on the free surface shape of oscillatory residual acceleration parallel to 
the axis of an isothermal liquid bridge have been made by Langbein [55]. Not suprisingly, his 
calculations suggest that the sensitivity of the surface is highest for disturbances with frequencies 
close to the natural frequency of the bridge, typically in the 0.001-10 Hz range. For the example 
given in [55], the maximum acceleration magnitudes that can be tolerated by a bridge less than 90% 
of the maximum stable bridge length (given by the Rayleigh limit) range from 5xl0‘ 8 to 5xl0' 3 g. 
The criterion used to determine the tolerable residual acceleration was that the amplitudes of all 
surface shape excitations with frequencies other than the resonance frequency be kept to 0.1% of 
the radius of the column. Examples of the results of this analysis are shown in Fig. 12. The curves 
represent the maximum tolerable axial acceleration for a given aspect ratio (Fig. 12a), and for a 
surface response with a given number of nodes (Fig. 12b). The results express the maximum 
tolerable frequency a 0 in terms of the g-jitter frequency to, the resonant frequency, co m , of the 
liquid zone associated with the m^ 1 node of the deformation of the zone surface, the viscosity of 
the liquid v and the spatial wavenumber 

a 0 < min[V (co 2 -ct&) - covq m ] , (9) 

where com = Y 9 m[( 9m R) 2 - l]/pR , and q m = (m+l)7t/H, R is the column radius and H the 
column height. Note that the global minimum, and local minima of the tolerable acceleration curve 
occur at disturbance frequencies equal to the natural frequencies of the system. 

Zhang and Alexander [91] have also analyzed this problem using a slender-body 
approximation. The problem of determining the axisymmetric response of the shape of the free 
surface of a cylindrical liquid column bounded by two solid regions is modelled by a 1-D system 
of non-linear equations. It is found that the sensitivity of the zone to breakage and shape changes 
depends on the static Bond number, aspect ratio and viscosity as well as the amplitude and 
frequency of the disturbance. The general trend is an increase in tolerable residual acceleration 
with increasing frequency. At the eigenfrequencies of the zone, however, there are strong 
deviations from this trend. At these frequencies this model also shows that the tolerable residual 
gravity level can be two orders of magnitude lower. These most sensitive frequencies have been 
found in the 10 _2 -10 _1 Hz range (see for example Fig. 13). For the cases examined, maximum 
tolerable residual gravity levels as low as 10' 6 g have been calculated. For higher viscosities the 
tolerable acceleration level is increased for all frequencies. The equilibrium shape, as determined by 
the steady background acceleration, has a prounounced effect on the response at low frequencies. 
A change in slenderness of the bridge markedly changes the sensitivity to residual acceleration as 
the Rayleigh limit is approached. It is clear that for liquid bridges, small magnitude accelerations 
cannot be ignored. 

For non-isothermal zones, the dimensionless groups introduced earlier in section 3.1 govern 

the problem. For molten silicon Sekerka and Coriell [89] have calculated the values of these 
groups. They are presented in Table 4. The relative magnitudes of the terms suggest that on earth 


12 



the zone shape will be dominated by the capillary effect (large Bo) for typical zone sizes, while 
under low-gravity conditions the zone shape will be influenced mainly by the temperature gradient 
along the interface. In addition (as discussed earlier) dynamic distortions may also be important. 
Considerations of the float zone crystal growth process introduce further complications which 
ultimately require independent analysis. In the previous cases, the surfaces of the inert solids 
between which the zone is suspended are usually planar. This is generally not the case for the 
crystal growth process. The actual shape depends primarily on the thermal conditions. The 
influence of the meniscus angle at either of the crystal-melt-vapor tri-junctions can also be 
important and has been analyzed by Heywang [93] and Coriell et al. [94]. In addition to zone 
failure due to the classical capillary instability (or Rayleigh instablity), a "dewetting” instability is 
found to be important. The instability arises when the angle between the meniscus and the melt- 
crystal interface is less than the wetting angle for the liquid melt on the solid crystal. For Si zones 
of greater than 1 cm radius, zone failure on earth occurs due to the Heywang instability [94,95]. 
Zone failure by this mechanism is unlikely under low-gravity conditions. 

Feuerbacher et al. [56] developed an ad hoc model to describe experiment sensitivity to 
residual accelerations. The model is based on the motion of a body in response to an applied 
sinusoidal oscillation. The governing equation is that of a damped harmonic oscillator. The 
tolerable residual accelerations are (as expected) found to increase with increasing frequency of the 
disturbance. Furthermore, there is a minimum tolerable acceleration at disturbance frequencies 
corresponding to the natural frequency of the system. 

The effect of time dependent accelerations on experiments involving drops and bubbles has 
been examined for the cases of sinusoidal and Heaviside-function disturbances [57]. The criterion 
used to define the sensitivity is that the displacement of the drop not exceed 10% of the drop 
radius. For a given set of physical properties, the admissible accelerations can be expressed as a 
function of frequency (Fig. 14). The tolerable acceleration levels exhibit trends similar to those 
depicted in Figs. 4-6. 

The influence of the gravity gradient, atmospheric drag and spacecraft attitude motion on the 
motion of a spherical solid particle immersed in viscous fluid has also been considered [5,6]. The 
extent and nature of the influence of the gravity gradient and the atmospheric drag on the trajectory 
of the particle are, for a given value of the Stokes drag coefficient, shown to be dependent on the 
characteristic orbital attitude, the magnitude of the atmospheric drag accelerations and the distance 
of the particle from the mass center of the spacecraft. 

Amin [59] has investigated the influence of g-jitter on heat transfer from an isothermal sphere 
maintained at a constant temperature greater than the ambient fluid temperature. The body force 
was taken to have a single frequency sinusoidal time dependence and a small amplitude. A 
significant result of the analysis is that the Reynolds stresses associated with the fluctuating flow 
result in a steady streaming motion. This has implications for both heat and mass transport in 
oscillating flows. For the low Pr fluids examined, it was concluded that buoyancy induced 
convection caused by high frequency g-jitter (with amplitudes 10* 2 g) will result in heat transfer 
conditions significantly different from pure conduction. However, it was found that in high Prandtl 
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number fluids with small kinematic viscosities, low frequency g-jitter will influence heat transfer. 

Analyses of heat transfer in (laminar) oscillating flows in cylinders and between parallel 
plates [27-33] have shown that at high frequencies the effective thermal diffusivity is 


Keff 06 


(0 Ax 2 
W 0 ’ 


while for low frequencies 


Keff «= co Ax 2 W D . 


( 10 ) 
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Here <o is the circular frequency, Ax is the cross-stream average displacement of a fluid element 
over half the period of the oscillation and W Q = L(o)/v)K 2 is the Wormsley number which 
represents the ratio of the viscous response time to the period of the oscillation. For a given 
frequency the heat transfer process was found to be most efficient for Pr = n /W 0 2 i.e. when the 
characteristic time, x k = L 2 /k associated with heat diffusion is equal to half the period of the 
oscillation. Fig. 15 illustrates this result for fluids with different diffusivities and L = 1 cm. 

2.3.3. Linear stability and gravity modulations 

The stability of a heated fluid layer of infinite lateral extent is affected by gravity modulation 
and has received a limited amount of attention for the case of sinusoidally modulated gravity [74- 
77]. The static case, where the gravity vector is perpendicular to the layer and parallel to the 
temperature gradient, is stable unless the critical Rayleigh number is exceeded. The effect of 
gravity modulation is to introduce stability for given combinations of modulation amplitude and 
frequency for situations which in the absence of modulations would exhibit instability. In contrast, 
for the case where the gravity vector and the temperature gradient are anti-parallel (which is stable 
for a steady gravity vector,) instability can occur due to gravity modulation. 

The work of Biringen and Peltier [77] focused on the behavior of a 3-D Rayleigh Benard 
system subject to spatial and temporal gravity modulation. They used a pseudo-spectral method to 
examine the system numerically and studied the effects of sinusoidal and random gravity 
modulations with the gravity vector oriented perpendicular to the fluid layer. It was found that 
when the base acceleration was reduced from lg to zero, the system parameters that had resulted in 
a synchronous response at 1 g led to a conductive state for the modulated zero g case. The 
spatially random modulations about a zero g base acceleration led to excitations of the 3-D 
temperature field. These excitations were less pronounced than for the 1 g case. For temporally 
random modulations about zero g, the cases considered that where excited at 1 g no longer 
exhibited excitation. 

The effect of oscillating accelerations and impulse forcing of an interface separating two 
immiscible fluids has recently been considered [78,79]. Both these analyses are restricted to time- 
dependent accelerations with zero mean which are oriented perpendicular to the interface between 
the fluids. The linear stability of a basic state characterized by a planar interface and an oscillating 
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pressure was examined. For the case of an arbitrary oscillating disturbance it is shown that the 
linear evolution of an infinitesimal perturbation of the interface is governed by a single integro- 
differential equation. Except in the limit of zero surface tension the effect of viscosity is shown to 
be small. As surface tension becomes dominant, resonant instabilities are excited at lower and 
lower wavenumbers. The interface is found to be more unstable to impulse type forcing than to 
sinusoidal forcing [79]. 

2.3.4. Direct numerical simulation 

The effects of low-gravity on the transport of heat and momentum have been examined in a 
number of articles [37, 60-71]. Robertson et al. [60,61] found that for convection in circular 
cylinders with azimuthal variations in temperature at the boundary and the gravity vector applied 
perpendicular to the cylinder axis, the intensity of convection follows the prediction of 
Weinbaum's first order theory [96] for low Rayleigh numbers. Weinbaum’s first order theory 
predicts a simple sinusoidal dependence of the maximum velocity as a function of orientation of the 
gravity vector. The effects of a variety of acceleration vectors (sinusoidal, cycloidal and linear 
periodic) on motion in three fluids corresponding to mercury, helium and water (10* 2 < Pr < 10) 
were studied by Spradley et al. [62]. They found that the steady mean part of the applied 
disturbance is more important than the oscillatory part (frequency =1 Hz) in determining the flow 
field and heat transfer. Kamotani et al. [37] solved a linearized approximation of the Boussinesq 
equations and investigated the effect of an applied acceleration consisting of a time mean part and 
an oscillatory part on the temperature and flow fields in a rectangular enclosure. They found that 
the thermal convection was predominantly oscillatory in nature. In addition they found that the 
acceleration perpendicular to the temperature gradient is the most important for the generation of 
fluid motion. 

The response of natural convection of a Boussinesq fluid in a cylinder to residual acceleration 
has been examined by Heiss et al. [64], and Schneider and Straub [97]. In the work of Heiss et 
al. the calculations are three-dimensional and involve gravity pulses and a rotating gravity vector. 
The ends of the cylinder are held at different temperatures. The cylinder walls are taken to be either 
perfectly insulating or perfectly conducting. The results are presented in terms of the relevant 
dimensionless parameters in the original papers. Figure 16 illustrates their results for impulse 
disturbances. The physical conditions were taken to correspond to: a cylinder of aspect ratio 1, 
with diameter = 1 cm, v = 5xl0" 3 cm 2 s" 1 , k= 7.14x10 -3 cm 2 s' 1 , (3 = 2.5xl0 -4 , and a 
temperature difference, AT, of 25° K. The magnitude of the acceleration pulse (duration =0.7 
seconds) corresponding to a Rayleigh number, Ra*, of 2000 is ~ 3xlO" 3 g and results in a 
maximum velocity of approximately 5.8 cm s' 1 . The decay time for the flow field is about 5.6 
seconds. Notice that for all the results, the decay time of the pulse is independent of the magnitude 
of Ra*. For a rotating gravity vector of a given magnitude, the maximum velocity of the system is 
independent of the dimensionless frequency f = f *d 2 /k, up to / ~ 10 whereupon it decreases 
rapidly with increasing frequency. For the fluid with the physical properties described in the 
previous example, / * corresponds to 7.14xl0' 2 s _1 . 

Schneider and Straub [97] examined the effect of Prandtl number on the response to g-jitter 
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for the same system as Heiss et al. They examined the response of the system to pulses of infinite 
length and sinusoidally varying pulses. For sinusoidal pulses they found for systems with Pr 
=0.71 and 7.0 (air and water) that the maximum velocity decreased with increasing frequency for 
dimensionless frequencies, F=/L 2 /k, greater than 10. At Pr=7 a maximum in the velocity was 
observed between F = 10 and 100 for the maximum amplitude pulse examined (Ra = 2xl0 5 , or I 
g I » 3xl(H g for a fluid with the properties of water, a 1 cm diameter cylinder and a temperature 
difference of 10 K). They determined that the system was least sensitive when the transient 
accelerations acted along the cylinder axis. 

The extent to which gravity causes buoyancy-driven fluid motion (and thus, solute 
redistribution) during directional solidification has been examined using numerical models of 
buoyancy-driven convection in cylindrical and rectangular geometries[9, 65-69, 84, 98-101]. Except 
for [9,65-69] these studies are restricted to axisymmetric situations in which a steady gravity vector 
is oriented parallel to the axis of a cylindrical ampoule. McFadden and Coriell [66] have 
undertaken 2-D calculations of the effects of time-dependent accelerations on lateral compositional 
variations during directional solidification. The gravitational acceleration was assumed to have a 
uniform magnitude and rotation rate. The amount of compositional nonuniformity was found to 
increase with decreasing rate of rotation (see Fig. 17). 

Two- and three- dimensional models of directional solidification from dilute gallium-doped 
germanium melts have been used to determine the sensitivity of crystal compositional uniformity to 
both time-dependent and steady residual accelerations [9,67,69]. The specific boundary 
conditions, thermo-physical properties of the melt, growth rates and ampoule size are all found to 
play a role in the determination of the experiment sensitivity. For a given set of operating 
conditions, it is found that at growth rates on the order of 6.5 mm s' 1 the orientation of the 
experiment with respect to the steady component of the residual gravity is a crucial factor in 
determining the suitability of the spacecraft as a means to suppress or eliminate unwanted effects 
caused by buoyant fluid motion. The worst case appears to be when the acceleration vector is 
parallel to the crystal interface. At growth rates on the order of microns per second, this orientation 
leads to compositional non-uniformities of 10-20% when the magnitude of the acceleration is of the 
order 10' 6 g. If, however, the growth rate is lowered by an order of magnitude, the nonuniformity 
is reduced significantly (down to 1- 5% in this case). Table 5 summarizes the results obtained for 
steady accelerations. 

The directional solidification (Bridgman-Stockbarger) process is also extremely sensitive to 
transient disturbances. The response of the system to a variety of impulses will be discussed. For 
example, a 3xl0* 3 g impulse of one second duration acting parallel to the interface of a growing 
crystal produces a response in the solute field which lasts for nearly 2000 seconds. Consequently 
lateral and longitudinal compositional variations occur over a length of nearly 6 mm in the grown 
crystal. 

The response of the solute field and the lateral nonuniformity to oscillatory accelerations 
varies from no response at all (at frequencies above 1 Hz with amplitudes below 10" 3 g) to a 
significant response at 10" 3 Hz at amplitudes on the order of 10~ 6 g. In addition, additive effects 
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were observed for combinations of a steady and a low frequency residual acceleration component. 
These additive effects gave rise to significant lateral and longitudinal non-uniformities in 
concentration. A number of different types of periodic disturbances were employed. Single 
frequency disturbances of the form g(t) = go + gn cos(2jt/ n t) were examined with g 0 = 0, 
^2(10)" 6 and v2(10)" 3 g, oriented parallel, perpendicular and at 45° to the crystal-melt interface. 
The frequency range examined was f n = 10 " 4 , 10' 3 , 10' 2 , 10' 1 , 1 and 10 Hz. For frequencies 
greater than 10* 2 Hz, there were no discernible effects on the solute fields. The velocity field did, 
however, respond to the oscillatory disturbances. For the case of 10" 3 Hz (at 5xl0' 6 g) the 
response of the solute field was significant. Lateral and longitudinal nonuniformity levels in excess 
of 15% were calculated. Figures 18-20 show the lateral nonuniformity as a function of time and 
highlights the additive effect of oscillatory and steady components of the residual acceleration. 

The effect of a multiple frequency disturbance is illustrated in Figure 20. The acceleration 
consists of steady and periodic contributions with the form: g(t) = go + gi cos(27tl0' 3 t) + g 2 
cos(2jtl0 -2 t ). Here I g 0 I = V2xl0'6 g, | gl | = 3 a/2x10-6 g and i g 2 I = 3^2x1 0'5 g. The 
magnitude of the compositional nonuniformity £, is seen to vary with the frequencies of the 
acceleration. More recently, it has been shown [67] that, particularly at high frequencies, it is 
important to consider the long time behavior of the system, because for high frequency 
disturbances ( « 0.1 - 1 Hz) the transient behavior of the system is more sensitive than the long 
time behavior. For example, the response to a 1 Hz 10' 2 g acceleration was found to exhibit a 
concentration nonuniformity on the order of 10% during the first 80 seconds whereupon it decayed 
slowly to 1% over a period of 2500 seconds. In addition, calculations have been made using 
sample acceleration time histories obtained on Spacelab 3 (SL-3). The results varied according to 
the type of SL-3 data that were input. For the example shown in Fig. 21 the response of the solute 
field was negligible, although the velocity field fluctuated on the order of the crystal growth rate. 
The system exhibited some response to transient disturbances (for example the structural response 
to thruster firings) that were superimposed on background low frequency acceleration. Low 
frequency (< 10' 1 Hz) acceleration components with amplitudes of 10" 3 g or more were also 
found to give rise to undesirable nonuniformity. 

Griffin and Mokatef [72,73] modelled the nature of melt convection in a Bridgman system 
subject to steady and unsteady axial acceleration. They also found that at conditions corresponding 
to low-gravity the response time of the melt velocity is proportional to the momentum diffusion 
time (i.e. it is controlled by the characteristic system length scale and melt viscosity). In accord 
with the general trend exhibited by most physical systems discussed in this article, they found that 
the velocity response to sinusoidal g-jitter decreased as the inverse square of the momentum 
diffusion time. 

A general conclusion that can be drawn from all attempts to characterize gravity-driven 
convective effects on directional solidification from two component melts is that the maximum 

lateral solute nonuniformity (radial segregation for the axisymmetric cases) occurs near the 
transition from diffusion dominated to convection dominated growth conditions [52,84,9]; that is, 
when convective velocities are of the same order of magnitude as the diffusive velocities. The 
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conditions under which this "transition" takes place will depend on the specific nature of the forces 
driving convection. The orientation of the steady component of the gravity vector is crucial in 
determining the magnitude of the gravity vector at which this transition occurs. Thus, for a given 
set of operating conditions the orientation of the gravity vector determines the suitability of a low- 
gravity environment for directional solidification experiments. 

Motivated by the growth of triglycine sulfate (TGS) crystals on Spacelab 3 [11], Nadarajah et 
al. [102] have examined the acceleration sensitivity of crystal growth from solution. Thermal and 
solutal convection were included in their analyses and the crystal growth rate was chosen as the 
sensitivity parameter for the response to convective transport. Simulations were carried out for 
steady, impulsive and periodic accelerations in order to determine tolerable acceleration levels. 
Long-time simulations of the experiment were conducted with steady background accelerations of 
10' 6 and 10' 5 g. Impulsive and periodic disturbances of higher magnitudes were imposed at 
intermediate points. For steady accelerations the system was shown to have diffusion dominated 
transport at 10' 6 g but for 10' 5 g a transition to convection dominated transport occurs. Figure 22 
summarizes the results for oscillatory disturbances. For each case, the numbers correspond to the 
ratio g(3Ac//v*, which is an estimate of the the ratio of buoyancy to inertial forces. Here (3 is the 
solute expansion coefficent, Ac the characteristic concentration difference,/ is the frequency of the 
acceleration and v* is the calculated characteristic fluid velocity. Clearly, for this example, a 
simple estimate alone will not suffice to determine the system sensitivity. The system is relatively 
stable to impulsive and periodic disturbances unless their magnitudes are very large. 

2.4. Summary and discussion 

The work described in this paper examines the sensitivity of a variety of space experiments to 
residual accelerations. In all the cases discussed, the sensitivity is related to the dynamic response 
of a fluid. In some cases the sensitivity can be defined by the magnitude of the response of the 
velocity field. This response may involve motion of the fluid associated with internal density 
gradients, or the motion of a free liquid surface. For fluids with internal density gradients, the type 
of acceleration to which the experiment is sensitive will depend on whether buoyancy-driven 
convection must be small in comparison to other types of fluid motion (such as thermocapillary 
flow), or fluid motion must be suppressed or eliminated (such as in diffusion studies, or 
directional solidification experiments). In the latter case, the experiments (for example diffusion 
experiments and directional solidification experiments) are sensitive to steady and low frequency 
accelerations (< 10~ 2 Hz). For experiments such as the directional solidification of melts with two 
or more components, determination of the velocity response alone is insufficient to assess the 
sensitivity. The effect of the velocity on the composition and temperature field must be considered, 
particularly in the vicinity of the melt-crystal interface. As far as the response to transient 
disturbances is concerned the sensitivity is determined by both the magnitude and frequency of the 
acceleration and the characteristic momentum and solute diffusion times. 

The directional solidification of doped germanium crystals has been carried out on Skylab 
[17,18,103,104] and the Apollo-Soyuz Test Project (ASTP) [105,106]. The latter involved the 
directional solidification of seeded Ge melts. Melts were doped with gallium and some were doped 
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with 1% Si and .001% Sb [105]. The results of these experiments revealed strong asymmetric 
non-uniformities in the space grown crystals. Lateral variations were also observed in samples 
grown under terrestrial conditions but were much less pronounced [25]. It has been argued that the 
asymmetric redistribution of the dopant can be ascribed to "barometric diffusion" of the solutes due 
to the acceleration gradient in the melt arising from the rotational motion of the spacecraft 
[105,106]. The basis of the argument, however, appears to ignore the presence of gravity gradient 
and atmospheric drag effects and does not explicitly account for the spacecraft attitude motions. It 
has been demonstrated [9] that melt convection can occur in response to a low magnitude steady 
residual acceleration, and that whenever the acceleration vector is not aligned with the ampoule axis 
strong asymmetries in composition can occur. Since there is a non-linear dependence of 
compositional uniformity on convection, there is always the possibility that the observed non- 
uniformities in the terrestrially grown samples will be smaller than those observed in space grown 
samples. Thus, convection caused by residual acceleration cannot be discounted as the origin of the 
non-uniformities. 

Thermocapillary flow experiments require that buoyancy effects are small in comparison to 
capillary effects. The estimated range of accelerations to which they are sensitive lies between 
5xlO" 3 Hz and 2 Hz for the most sensitive cases. Unless the magnitude of the effective buoyancy 
force (expressed through the Grashof number) associated with a given acceleration is on the order 
of or greater than Rs(l+ St 2 ) 1 / 2 , buoyancy effects are not likely to be significant. It is possible that 
more subtle effects of the acceleration environment might occur (for example influencing the 
transition from steady to oscillatory flow) but this cannot be determined by order of magnitude 
analysis. It is more likely that thermocapillary flow experiments will be most sensitive to the 
motion of the free surface in response to g-jitter. 

The response of free surfaces to g-jitter has been examined through analytical and numerical 
models. For cylindrical liquid bridges this has been undertaken only for axial accelerations. The 
shape of the surface is most responsive to disturbances close to the natural frequency of the bridge 
and its harmonics. The sensitive frequencies lie in the 5x1 0' 3 -10 Hz range, with maximum 
tolerable magnitudes from 10" 8 - 10' 3 g depending on the sensitivity criterion and the particular 
model employed. 

Figure 23 depicts the range of measured acceleration amplitudes (courtesy of Dr. H. 
Hamacher, DLR, and M. J. B. Rogers, University of Alabama in Huntsville) as a function of 
frequency for several “shuttle” missions. It should be remembered that in addition to the 
accelerations shown here there are gravity gradient and atmospheric drag accelerations with 
magnitudes on the order of 10' 6 -10" 7 g which can be steady (for gravity gradient stabilized type 
attitudes) or have frequencies that are twice the orbital frequencies (for a solar inertial attitude). In 
Fig. 24 we depict these together with g-tolerance curves selected from the examples discussed 
earlier. There is some overlap between the sensitivity curves and measured accelerations. Even if 
the fact that most of the curves are obtained for order of magnitude estimates is taken into account, 
it is clear that a careful evaluation of the interaction between experiments and the residual 
acceleration environment is necessary. Once the experiment sensitivity has been determined 
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(which may sometimes require a more detailed analysis than any discussed in this review) it should 
be possible for the experimenter to choose operating conditions which will minimize the effect of 
the predicted acceleration environment. Since flight opportunities are limited, it is particularly 
important that the chances of success for any given experiment be greatly improved by detailed 
modelling prior to flight. This will allow investigators to better anticipate undesirable effects of the 
prevailing residual acceleration conditions, and will ultimately lead to a better overall understanding 
of the physical system. 
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2.6. Tables 

Table 1. Changes in residual acceleration components [|ig m*l] along the local vertical (Xi) and 
perpendicular to the oribital plane (X 2 ) associated with the gravity gradient stabilized 
attitude, and estimated atmospheric drag accelerations, a and b are, respectively, best and 
worst drag estimates for the “Shuttle” (1) and NASA’s Space Station (2). From [5] 


Altitude [km] 

Xi 

x 2 

Drag la 

Drag lb 

Drag 2a 

Drag 2b 

275 

0.414 

-0.137 

2.10 

0.360 


3.1 

300 

0.410 

-0.136 

1.20 

0.210 

6.1 

1.8 

325 

0.405 

-0.134 

0.69 

0.120 

3.6 

1.1 

350 

0.401 

-0.133 

0.41 

0.070 

2.1 

0.6 

400 

0.392 

-0.130 

0.24 

0.025 

0.7 

0.2 

450 

0.383 

-0.127 

0.05 

0.009 

0.4 

0.1 

500 

0.375 

-0.124 

0.02 

0.001 

0.8 

0.02 
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the composition at the interface. 


g-orientation 
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61 
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1.25 
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25 
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21 
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145 
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0.27 

2.27 
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1(b) 
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[92]. The length scale is taken to be the zone radius R; Pt is estimated at 1.6 x 10~ 4 K' 1 . For 
space, gs is taken as ICHg. AT is assumed to be 10 C°. (After Sekerka and Coriell [89]). 
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Table 5. 


Compositional nonuniformity % [%] for Ge:Ga [9,80]; \ - (c smax -c sm i n )/c sa v, 
where c s denotes the crystal composition at the melt-crystal interface. 


Residual Orientation 


Ampoule Width [cm] 


Acceleration 


1 0.5 

2.0 

Magnitude 


Growth Rate [mm s* 1 ] 


[g] N e g 

6.5 

0.65 6.5 

6.5 

A) 10- 4 

80 



10' 5 

92.7 

11.9 12.0 



70.9 

11.3 



6.4 

0.95 


5 ( 10)'6 

3.2 
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54.2 



io - 6 

11.3 

2.0 



8.0 




0.7 

0.0 


B ) IO ' 5 

22.6 


64.5 

10-6 

2.3 



eg is the unit vector parallel to g, N is the normal vector to the crystal. L = 1 cm for all cases. Values in 
parentheses indicate 3-D results. All calcualtions were undertaken using the Thermo-physical properties of gallium- 

doped germanium [80] and A) and B) refer to the operating conditions listed below. 


Operating conditions 




A ) Hot zone temperature (Th) 


1331 K 


Distance between inlet and interface (L) 


1.0 cm 


Height of adiabatic zone (Lg) 


2.5 mm 


Ampoule width (diameter) 


1.0 cm 


B) Hot zone temperature (Th) 


1251 K 







2.7 . Figures 


Fig 1. 

Fig 2. 
Fig. 3 

Fig. 4 

Fig. 5 
Fig. 6 
Fig. 7 

Fig. 8 

Fig. 9 

Fig. 10 

Fig. 11 

Fig. 12 


An example of a relatively “quiet period” measured on SL-3 (< 4xl0 _3 g, g = 9.8 m 
s' 2 ), a) Total acceleration array, b) Combined amplitude spectrum, note 
components in the 5-6 Hz and 7-8 Hz, 10-12 Hz, 17Hz and 30-35 Hz ranges. 
After Rogers and Alexander [47]. 

“Window” containing three thruster pulses, a) Total acceleration array, maximum 

magnitude ~lxl0' 2 g. b) Combined amplitude spectrum, dominant frequency 
components 1 1 Hz and 17 Hz. After Rogers and Alexander [47]. 

Estimated tolerable residual accelerations (in units of g = 9.8 m s' 2 ) as a function 
of frequency for: 1) a fluid physics experiment involving a temperature gradient, 2) 
a crystal growth experiment and 3) a thermo-diffusion experiment. After Monti et 
al. [38]. 

Estimated tolerable residual accelerations (in units of g = 9.8 m s" 2 ) for 
semiconductor and metal solidification experiments. The sensitivity parameter, is 

longitudinal segregation, Ac/c . Pe s = 5, temperature gradient = 50K cm -1 , k=0.1. 
After Monti et al. [36]. 

Estimated tolerable residual acceleration (in units of g = 9.8 m s' 2 ) ranges for 
diffusion experiments. After Monti et al. [36]. 

Tolerable acceleration limits for a thermocapillary flow experiment expressed in 
terms of Grashof (Gr), Strouhal (St) and Surface Reynolds (Rs) numbers. 
Tolerable acceleration limits (in units of g = 9.8 m s' 2 ) as a function of frequency 
for the thermo-capillary experiments with properties listed in Table 2. a) Gravity 
forces 1% of thermo-capillary forces, b) gravity forces 10% of thermo-capillary 
forces. 

A plot of GrSc vs. Pe g for the results of numerical simulations for steady 
accelerations parallel to the crystal-melt interface [80]. The solid and dashed lines 
delineate the following transport regimes: I a - c convective regime; II advective- 
diffusive regime [83]. Note that the adiabatic zone length has been used to calculate 
the Grashof and growth Peclet numbers, and that x 0 =0.27. Numbers in 

parentheses indicate the value of compositional nonuniformity % = lc sma x - 

Csmin 1/csavX 100%, where c s denotes the crystal composition at the melt-crystal 
interface, associated with each of the numerical results. 

Transient convection in large and small test cells due to residual accelerations 
caused by astronaut motion after Dressier [54]. Here x is the angular displacement, 

v the velocity and ti /2 is the decay time given by 0.68xR 2 /v. 

Stability limits of an isothermal liquid bridge of volume V and length L suspended 
between discs of radius R and subject to a steady axial acceleration defined by a 

Bond number B 0) after Martinez et al. [41]. The angle ()> is defined in the inset. 

The maximum stable zone length of an isothermal liquid zone suspended between 
inert solids of radius R as a fucntion of Bond number: Comparison of theory and 
experiment. After Coriell et al. [90]. 

Tolerable residual accleration (in units of g = 9.8 m s' 2 ) for experiments 
investigating resonance frequencies of liquid columns subject to axial acceleration , 

following Langbein [55]. a) sensitivity for different aspect ratios A = L/2 tiR, b) 
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Fie 13 


Fig. 14 


Fig. 15 
Fig. 16 


Fig. 17 


Fig. 18 


Fig. 19 


Fig. 20 


Fig. 21 


change in sensitivity with number of nodes in the deformation for A = 0.9999. 
Curves of tolerable acceleration (in units of g = 9.8 m s* 2 ) versus frequency for a 
maximum liquid bridge shape change of 10% at B o =0.002, g 0 =1.42xl0‘ 5 g and 
C=0.001. The solid curve, the dotted curve and the dotted-dashed curve are the 
results for aspect ratios L/2R=2.6, 2.826 and 3.024, respectively. From Zhang 
and Alexander [91]. 

Threshold intensity g* = g T R 4 /v 2 d x as a function of normalized frequency S = 

coR 2 /2v after Dewandre [57]. Here g T is the maximum tolerable acceleration, d x the 
admissible displacement and r is the ratio of the drop density to that of the 
surrounding fluid. 

Plot of effective diffusivity, Keff vs. Wormsely number W 0 , after Kurzweg [31]. 

Dimensionless maximum velocity caused by impulse disturbances as a function of 
Fourier number F 0 = tK 2 /d, after Heiss et al. [64]. Here t is the time [s], k is the 

thermal diffusivity [cm 2 s* 1 ], d is the diameter of the cylinder [cm]. The parameter 
Ra* is the Rayleigh number associated with the disturbance which is of 

dimensionless duration F 0 *= 5.5(10)' 3 . Ra* = Pt gL 3 AT/vtc, where L is the 

cylinder length [cm], AT is the temperature difference between the ends of the 

cylinder [K], px is the coefficient of thermal expansion [K _1 ], g is the magnitude 

of the acceleration [cm s -2 ] and v is the kinematic viscosity [cm 2 s" 1 ]. 

Maximum lateral solute segregation Ac/c«> (where Coo is the far-field solute 
concentration in the melt) in the crystal as a function of the period of the 
gravitational oscillation, after McFadden and Coriell [66]. 

Lateral nonuniformity in composition, £,, plotted as a function of time for an 
oscillatory residual acceleration with a maximum magnitude of 3^2(1 0)' 6 g and a 
frequency of 10' Hz, acting parallel to the crystal-melt interface [67]. The initial 
state was purely diffusive. 

Lateral nonuniformity in composition, plotted as a function of time for a residual 
acceleration consisting of a steady part with magnitude V 2(10)' 6 g and an 
oscillatory part with a maximum magnitude of 3^2(1 0)' 6 g and a frequency of 10' 3 
Hz, acting parallel to the crystal-melt interface [67]. The calculation was started 
from a steady flow associated with a a/ 2(10) -6 g acceleration acting parallel to the 
interface. 

Lateral nonuniformity in composition, plotted as a function of time for a multi- 
component disturbance consisting of a steady low g background plus two periodic 
components: g(t) = g 0 + gi cos(27tlO"3t) + g 2 cos(27tl0~2t), where I g 0 I = 
V2(10)'6, | gi | = 3V2(10)'6 and|g 2 I = 3V2(10)'5 [67], The calculation was 

started from a steady flow associated with a V2(10)" 6 g acceleration acting 
perpendicular to the interface. 

Lateral nonuniformity % I as a function of time for accelerations taken from a 
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sample time series (see inset) constructed from data obtained on Spacelab 3 [69]. 
The acceleration consists of a repeated “noise” segment (length 10 seconds, 10 _1 
Hz </< 1 Hz) and a thruster firing event (length 10 seconds, 10' 1 Hz </< 13 
Hz). The latter is introduced at 10 and 80 seconds. 

Fig. 22 Sensitivity of crystal growth from solution to various periodic disturbances. The 
solid line separates tolerable responses (10% or less growth rate fluctuations) from 

intolerable ones. For each case, the numbers correspond to the ratio g[3Ac//v* , 
which is an estimate of the the ratio of bouyancy and inertial forces. After 
Nadarajah et al. [102]. 

Fig. 23 Residual accelerations measured on orbit (in units of g = 9.8 m s' 2 ). (Data 
courtesy of Dr. H. Hamacher, DLR, and M.J.B. Rogers Univ. of Alabama in 
Huntsville.) The Spacelab 3 (SL-3) data are restricted to measured frequencies 
between 10* 3 and 1 Hz. The points corresponding to Sled, Hop, Drop, Treadmill, 
Quiet, FPM ops, Stowage, refer to activities and experiments on the D-l mission. 
FPM ops stands for Fluid Physics Module operations, P RCS and V RCS refer to 
the primary and vernier thrusters. The drag and rotation entries correspond to 
accelerations arising from slow variations in atmospheric drag during an orbit and 
attitude changes involving rotation. 

Fig. 24 Acceleration sensitivity curves from selected examples discussed earlier 
superimposed on Fig. 23. Curves 1 and 3 correspond to the cases n=l, A=0.9999 

and n=l, A =0.9 in Fig. 12 a, curve 2 corresponds to n=2 and A=0.9999 in Fig. 
12 b. Curve 4 corresponds to the most sensitive case of zone shape change in Fig. 
13, and curves 5 and 6 respectively correspond to the semi-conductor and metal 
melt growth experiments from Fig. 4. Curve 7 represents the thermo-diffusion 
experiment of Fig. 5. 
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Amplitude Spectrum [1 0" 4 g] Residual Acce 






Log Tolerable Acceleration [g] 



Fig s 4 




Self-diffusion 


Inter-diffusion 


Thermo-diffusion 






Grashof Number (Gr) 



Strouhal Number (St) 


Fig.6 



Log Tolerable Acceleration [g] 



Log Frequency [Hz] 


Fig. 7 a 



Log Tolerable Acceleration [g] 




log GrSc 



log V g L a /D = log Pe g 


Fsg. 8 



1.2x10' 3 cm s' 1 


x= 0.06 



D = 10 cm 


* 1/2 170 8 


Time (sec) 



= 1 . 2 x 10 3 cm s 1 


kx = 0 . 28 ' 

'K 

\ 

\ 

\ 

\ 


Small 


D = 2 cm 


*i/i 7 8 



Fig, 9 



o 





Log Frequency [Hz] 


Fig. 12 a 








Log g 






n 

ui' 


Fourier Number 


0.04 




• Period 


Fig. 17 







26 



9 - 19 





Fig. 20 





Acceleration Amplitude [g] 





Log Acceleration [g] 


TS 

s 

m 

JO 

w 



Log Frequency [Hz] 


□ Sled 

* Hop 

• Drop 

° Treadmill 

* Quiet 

0 FPM ops 

+ Stowage 
A PRCS 

° V RCS 

O Drag 

v Rotation 

• SL 3 < 1 Hz 



Log Acceleration [g] 


“n 

BUSES £2 

<Q 

85 

to 

& 



Log Frequency [Hz] 



Page intentionally left blank 


3. A NUMERICAL ANALYSIS OF LOW GRAVITY 

TOLERANCE OF THE BRIDGM AN-STOCKB ARGER 

TECHNIQUE : TIME-DEPENDENT ACCELERATION 

3.1. Introduction 

In this work we extend our investigation of the sensitivity of directional solidification by 
the Bridgman-Stockbarger technique to residual acceleration and examine the effects of periodic 
accelerations and acceleration profiles synthesized from data recorded on Spacelab 3 (SL-3). The 
effects of low-gravity on the transport of heat and momentum have been examined in a number of 
articles [1-13]. Robertson et al. [2,3] found that for convection in circular cylinders with azimuthal 
variations in temperature at the boundary and the gravity vector applied perpendicular to the 
cylinder axis, the intensity of convection follows the prediction of Weinbaum's first order theory 
[14] for low Rayleigh numbers. Weinbaum's first order theory predicts a simple sinusoidal 
dependence of the maximum velocity as a function of orientation of the gravity vector. The effects 
of a variety of acceleration vectors (sinusoidal, cycloidal and linear periodic) on motion in three 
fluids corresponding to mercury, helium and water (10 -2 < Pr < 10) were studied by Spradley et 
al. [4]. They found that the steady mean part of the applied disturbance is more important than the 
oscillatory part (frequency =1 Hz) in determining the flow field and heat transfer. Kamotani et al. 
[1] solved a linearized approximation of the Boussinesq equations and investigated the effect of an 
applied acceleration consisting of a time mean part and an oscillatory part on the temperature and 
flow fields in a rectangular enclosure. They found that the thermal convection was predominantly 
oscillatory in nature. In addition they found that the acceleration perpendicular to the temperature 
gradient is the most important for the generation of fluid motion. 

The response of natural convection of a Boussinesq fluid in a cylinder to residual acceleration 
has been examined by Heiss et al. [6], and Schneider and Straub [7]. In the work of Heiss et al. 
the calculations are three-dimensional and involve gravity pulses and a rotating gravity vector. The 
ends of the cylinder are held at different temperatures. The cylinder walls were taken to be either 
perfectly insulating or perfectly conducting. Schneider and Straub [7] examined the effect of 
Prandtl number on the response to g-jitter for the same system as Heiss et al. They examined the 
response of the system to pulses of infinite length and sinusoidally varying pulses. For sinusoidal 
pulses they found for systems with Pr =0.71 and 7.0 (air and water) that the maximum velocity 
decreased with increasing frequency for dimensionless frequencies, F=/L 2 /k, greater than 10. At 
Pr=7 a maximum in the velocity was observed between F = 10 and 100 for the maxium amplitude 
pulse examined (Ra = 2xl0 5 , or I g ! ~ 3xl0' 4 g for a fluid with the properties of water, a 1 cm 
diameter cylinder and a temperature difference of 10 K). They determined that the system was least 
sensitive when the transient accelerations acted along the cylinder axis. 

The extent to which gravity causes buoyancy-driven fluid motion (and thus, solute 
redistribution) during directional solidification has been examined using numerical models of 
buoyancy-driven convection in cylindrical and rectangular geometries [7-11,15-20]. Except for [ti- 
ll, 20] these studies are restricted to axisymmetric situations in which a steady gravity vector is 
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oriented parallel to the axis of a cylindrical ampoule. McFadden and Coriell [9] have undertaken 2- 
D calculations of the effects of time-dependent accelerations on lateral compositional variations 
during directional solidification. The gravitational acceleration was assumed to have a uniform 
magnitude and rotation rate. The amount of compositional nonuniformity was found to increase 
with decreasing rate of rotation (see Fig. 17). 

Griffin and Mokatef [21,22] modelled the nature of melt convection in a Bridgman system 
subject to steady and unsteady axial acceleration. They also found that at conditions corresponding 
to low-gravity, the response time of the melt velocity is proportional to the momentum diffusion 
time (i.e. it is controlled by the characteristic system length scale and melt viscosity). They found 
that the velocity response to sinusoidal g-jitter decreased as the inverse square of the momentum 
diffusion time. 

Two- and three- dimensional models of directional solidification from dilute gallium-doped 
germanium melts have been used to determine the sensitivity of crystal compositional uniformity to 
both time-dependent and steady residual accelerations [10,11,20] (see also section 5.1). The 
specific boundary conditions, thermo-physical properties of the melt, growth rates and ampoule 
size are all found to play a role in the determination of the experiment sensitivity. For a given set 
of operating conditions, it is found that at growth rates on the order of 6.5 mm s" 1 the orientation 
of the experiment with respect to the steady component of the residual gravity is a crucial factor in 
determining the suitability of the spacecraft as a means to suppress or eliminate unwanted effects 
caused by buoyant fluid motion. The worst case appears to be when the acceleration vector is 
parallel to the crystal interface. At growth rates on the order of microns per second, this orientation 
leads to compositional non-uniformities of 10-20% when the magnitude of the acceleration is of the 
order 10‘ 6 g. If, however, the growth rate is lowered by an order of magnitude, the nonuniformity 
is reduced significantly (down to 1- 5% in this case). Table 5 in section 5.1 summarizes the 
results obtained for steady accelerations. 

3.2. Formulation 

The basic model and governing equations were introduced in paper I. The physical 
situation under investigation involves an ampoule containing a dilute two-component melt. 
Directional solidification takes place as the ampoule is translated through fixed "hot " and "cold" 
zones. The zones are separated by a thermal barrier [23] which is modelled using adiabatic 
sidewalls. Translation of the ampoule is modelled by supplying a doped melt of dilute bulk 
composition Coo at a constant velocity Vg at the top of the computational space (inlet), and 
withdrawing a solid of composition c s = c s (x,t) from the bottom. The crystal-melt interface is 
located at a distance L from the inlet. The temperature at the interface is taken to be T m , the melting 
temperature of the crystal, while the upper boundary is held at a higher temperature Th. Since we 
wish to confine our attention to compositional non-uniformities caused by bouyancy convection, 
rather than variations resulting from non-planar crystal-melt interfaces, the interface is held flat in 
our model. 
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The governing equations are cast in dimensionless form using L, k/L (k is the melt's 
thermal diffusivity), p m K 2 /L 2 , Th-T m , and Coo to scale the lengths, velocity, pressure, temperature, 
and solute concentration. The dimensionless equations governing momentum, heat and solute 
transfer in the melt are then 


-=r- + (gradu)u = PrAu + GrPr 2 g, 
dt 

div u =0, 

~ + u • gradG = A0, 
dt 

+ u • grads) = As, 


( 1 ) 

( 2 ) 

(3) 

(4) 


where, u(x) represents the velocity, 0 = (T(x) - T m )/(Th - T m ) the temperature and s = c m (x)/c” 
represents the solute concentration. Pr, Ra and Sc are defined in section 5.1. The magnitude of g 
used in the Rayleigh number, Ra, will be the magnitude of the actual acceleration. For the 
calculations presented here we have used the thermo-physical properties corresponding to gallium- 
doped germanium melts (see section section 2.6 p. 28 and section 5.1). The term g in (1) 
has a specifies the orientation of the gravity vector. The Rayleigh number is taken to be the value 
of Ra at the earth’s surface, thus the magnitude of g represents the acceleration magnitude relative 
to 1 g. 

The following boundary conditions apply at the crystal-melt interface: 

+ u • grads) = As, (5) 

where N is the unit vector normal to the planar crystal melt interface. We define the measure of 
compositional nonuniformity in the crystal at the interface to be the lateral range in concentration 
given by 


^(t) = , C s max-Csmin %> (6) 

c sav 

where c s is the (dimensional) solute concentration in the crystal, and c a v is the average 
concentration. The following boundary conditions are applied at the "inlet" (y=0) 


0 = 1, u-N = 
At the side walls the conditions are 


Pe g Sc 

Pr 


<j, uaN = 0, i^ = ?^l(s-l). 
dy Pr 


(7) 


u-N 


Pe g Sc 

Pr 


0 ? 


u-e w = 0, grads-e w = 0, 


( 8 ) 
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with 9 = 1 in the isothermal zone and grad0-e w in the adiabatic zone. Here e w is the normal to 
the ampoule wall. 

In an actual experiment, owing to the finite length of the ampoule, there is a gradual 
decrease in length of the melt zone during growth. In this model transient effects related to this 
change are ignored. This assumption is referred to as the quasi-steady assumption and is frequently 
used in melt-growth modelling[15, 19,20]. It is thus assumed that the ampoule is sufficiently long 
for these effects to be negligible, a condition which should be checked from case to case. 

3.2.1. Solution Method 

The governing equations were recast in the stream-function vorticity formulation and 
solved using a pseudo-spectral Chebyshev technique [20,24] which incorporated the influence 
matrix method [25,26]. This allowed us to avoid the artificial compressibilty method used to 
examine steady and transient accelerations and, thus, avoid iterating to satisfy (2) at each time 
step. 


3.2.2. Types of transient acceleration examined 

The form of accelerations examined was taken to be 

N 

g(t) = X g n sin(27t/n t +<)> n )k, (1) 

n=l 

where/n is the frequency, k is a unit vector parallel to the crystal melt interface, <j) n is the phase 
and g n represents the amplitude of the n 1 * 1 acceleration component. We examine single, triple 
component accelerations as well as accelerations synthesized from accelerations recorded on SL-3. 

The computations were performed on a CRAY XMP computer. The results of the 
calculations were compared with results obtained using the primitive variable formulation and 
artificial compressibility method used in paper I. The results were found to be practically identical. 
For the time dependent cases examined we systematically varied the number of spatial points as 
well as the time step size. In all cases a 33x33 Chebyshev points were found to be sufficient. The 
results presented here were calcualted using 41x41 Chebyshev points. 

3.3. Results 

We examined single and multiple frequency acceleration and synthesized accelerations from 
sample time series obtained on SL-3. Our results showed that, qualitatively, the transient response 
of the system to single and multiple frequency accelerations is independent of the type of 
disturbance. If examined over short time periods, that is t < L 2 /D, the sensitivity can be a factor of 
ten higher than at later times. 

a) Single frequency oscillatory disturbances 
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Table Igives the short time response of % for single frequency oscillatory disturbances. We 
shall see later that results obtained for longer simulation times indicates that a better estimate of the 
long time sensitivity is obtained after times on the order of L 2 /D 

We have examined the system response to disturbances in the frequency range 10' 3 -1 Hz. 
The most sensitive response was found to be for the longer period disturbances. The development 
of the velocity field is characterized by a plot of u max and v max as a function of time in Fig. 1. Note 
that for the 10' 3 Hz disturbance the largest magnitudes are attained at 250, 750 1250 and 1750 
seconds. The velocity field is in phase with the residual acceleration. The short time sensitivity is 
greatest at the end of the velocity transient. On the other hand, the compositional nonuniformity 
(see Fig. 2) exhibits a local maximum at about 600, 1000, and 1500 seconds. For the first two 
periods of the 10- 3 Hz disturbance the solute field is not in phase with the velocity field. The 
solute field at 500, 1000, 1500 and 2000 seconds for this case is shown in Fig. 3. 

At higher frequencies the response of the velocity field changes character. The difference 
between the behavior for 1 Hz, 10' 1 Hz and 10‘ 2 Hz is marked. The maximum velocity exhibits a 
short transient of between 20-50 seconds before it begins to oscillate regularly. Note that, in 
comparison to the 10" 3 Hz case, the maximum value of % for the 10' 2 and 10' 1 Hz cases was 
reached between 200-300 seconds, whereupon it decays significantly. Figures 4-7 show the 
maximum velocities and concentrations as a function of time for the 10 -2 Hz and 10" 1 Hz cases. 
Note the appearance of two local maxima of ^ per period at about 1200 seconds for the 10 Hz 
case. 

Fig. 8 shows the velocity field at .25 s, 0.5s, 0.75 s and Is for a 10' 2 g 1 Hz acceleration 
oriented parallel to the crystal surface. Note that in all the velocity plots, the velocity vectors have 
been rescaled such that the velocity at the inlet and the crystal boundary is always of magnitude 
6.5x1 0‘ 4 cm s _1 . During first half of the cycle the acceleration vector points from right to left. The 
maximum magnitude occurs at 0.25 s. The hot melt moves anti-parallel to the acceleration vector. 
The velocity reaches a maximum at 0.25 seconds and then decreases as the acceleration magnitude 
decreases. At 0.75 seconds the acceleration vector has reached a value of 10' 2 g but is now 
oriented from left to right. The fluid motion at this time, however, still has the same sense as the 
motion during the first half of the cycle, although the magnitude is diminished. At 1 second, 
significant convection occurs only at the boundaries. The flow direction has reversed and a weak 
clockwise circulation has been established. This is short-lived since the acceleration direction then 
changes for the next 0.5 seconds. This behavior repeats itself for some time, although the relative 
strength of the clockwise rotation increases at the expense of the counter-clockwise flow. This is 
seen in Fig. 9. Here the component (u,v) maximum velocity magnitude are seen to fluctuate with 
twice the period of the driving force. During the first second, a strong counter-clockwise cell is 
established followed by a weak clockwise cell. This is indicated by the appearance of two velocity 
maxima, the first representing the counter-clockwise cell which is approximately 2.5 times the 
magnitude of the clockwise cell. During subsequent cycles, the difference between these two peaks 
decreases. At about 80 seconds the peaks have approximately the same magnitude. Fig. 10 shows 
the component maximum velocities from 802 to 808 seconds. The behavior is quite regular and 
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about 7t/2 out of phase with the driving force. Notice the slight difference between adjacent peak 
velocity magnitudes. If the two maxima per period are the same, then no net flow will have 
occurred in the system over that period. Careful examination of Fig. 10b reveals a barely 
perceptible difference in the two maxima.This reflects a small net circulation with velocities on the 
order of the growth velocity, 

Fig. 11 depicts the solute concentration field at 1, 6, 200 and 800 seconds. At 1 second the 
solute field appears to be slightly disturbed, at 6 seconds the isoconcentrates away from the crystal- 
melt interface are clearly distorted, although the nonuniformity at the surface is negligible. At 200 
seconds a distinct nonuniformity has been established. At 800 seconds the magnitude of the 
nonuniformity has decreased but is still significant. The variation of t, with time is depicted in Fig. 
12. Here, £ reaches a value of 19% at about 300 seconds and decays to almost 6 % after 3000 
seconds. As shown in the inset b) ^ is modulated slightly with a frequency of 1 Hz. 

b) multiple frequency oscillatory disturbance 

In order to assess the contribution to solute transport of a combination of oscillatory 
acceleration components, we examined several cases of a three component disturbance. As for the 
single component cases, for the range of amplitudes chosen, the greatest sensitivity was to the low 
frequency component of the disturbance. In these cases the lowest frequency was 10' 3 Hz. 

Fig. 13 depicts the velocity field at 3,6,8 and 10 seconds in response to a disturbance with 
frequency components 10' 2 , 10' 1 , and 1 Hz and amplitudes 10 -4 , 10‘ 3 and 10" 2 g. Development of 
the velocity field during the first 10 seconds is illustrated by Fig. 14. Fig 15 gives the 
corresponding solute concentration fields. For this case the maximum solute nonuniformity was 
almost 60% . While the effects of the individual components are additive, the higher frequency 
accelerations cause only relatively small modulations of the basic trend. For a disturbance with 
frequency components 10' 2 , 10' 1 , and 1 Hz and amplitudes 10' 5 10~ 4 and 10' 3 g, the maximum 
velocity and composition nonuniformity as a function of time are shown in Figs. 16-17. As for the 
case of the single frequency disturbances, the velocity field gives rise to a counter clock-wise net 
circulation during the transient phase. The maximum value of is 6%. This is attained at around 
250 seconds and decays to value which fluctuates between 0.18 and 0.65% after some 3000 
seconds have elapsed. Clearly, this case is only marginally sensitive for long times. In addition, it 
is interesting to note the behavior of the nonuniformity fluctuation which is sinusoidally modulated 
during its overall decay after 200 seconds. The most obvious modulations occur with a frequency 
of 10' 2 Hz (Fig. 18a) and these are not manifested until 80 seconds have elapsed. Fig. 18b shows 
the long time behavior of % from 1200 seconds onward. The mean value of ^(averaged over 100 
seconds) continues to decay until about 2000 seconds whereupon it fluctuates about a value of 0.4 
%. Notice also that at 2000 seconds an additional local maximum in £ appears such that the 
modulation period is 50 seconds. This merely reflects that the composition field is finally 
responding to the basic flow which had been established some 1900 seconds earlier. Now the flow 
results in two maximum values of £, per period. The second maximum, which is related to the 
clockwise motion of the melt, is suppressed until tp reaches a steady value. The velocity at a point 
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x=0.7, y=0.7 as a function of time from 800-1000 seconds is shown in Fig. 19. For this case, the 
10‘ 2 Hz component of the acceleration is responsible for the largest fluctuations. 

In all the cases examined, the velocity field exhibited a greater response than the solute field. 
For single frequency oscillatory disturbances, the degree to which the fluid motion affected the 
time-dependent redistribution of solute in the melt is dependent on the characteristic time for the 
fluid to reach an oscillatory state compatible with the driving force. If the the Schmidt number is 
large the solute field will be disturbed during the velocity transient. At the time the velocity 
transient ends and the fluid velocities oscillate between fixed values at the frequency of the driving 
force, the solute field will still be in a transient state. In fact, since during the transient there will 
have been net rotation of the fluid in one direction, the composition nonuniformity will reach a 
maximum toward the end of the velocity transient. Once the fluid behaves regularly the solute field 
slowly adjusts to this state. The composition nonuniformity gradually reduces and reaches a 
significantly smaller value after an elapsed time on the order of the characteristic diffusion time Tp 
(more than 2000 seconds for the Ge:Ga system examined here). 

c) Acceleration synthesized from SL-3 sample time series 

Of the cases examined all but one are on the margin of sensitivity for the model Ge:Ga 
system. The results are presented in Figs. 20-24. The compositional nonuniformity was insensitive 
(^<0.1 %) to the disturbance shown in Fig. 20, although the maximum velocity fluctuated on the 
order of the growth rate (see Fig. 21). For the acceleration depicted in Fig. 20, the amplitudes of 
the component accelerations are given in Table 2. In order to find the margin of sensitivity of a 
disturbance with these frequency components we increased the amplitudes by a factor of lO.Fig. 
22 shows the time variation of u ma x and v max is shown. Note that the sample series used to 
synthesize the disturbance has a period of 10 s which is evident in Figs.21, 22 and 23. As might 
be expected, the tenfold increase in the component amplitudes resulted in a tenfold increase in the 
maximum value of 

In Fig. 24 the response to a disturbance synthesized from spacelab-3 acceleration data is 
shown. The disturbance consists of two repeated sample time series (a) (0-10 seconds, 20-80 
seconds and from 90 seconds on, combined with sample time series (b) (from 10-20 and 80-90 
seconds. Table 3 gives the component frequencies and amplitudes of disturbance (a). The sample 
series (b) represents a thruster firing event. For the amplitude-frequency components present in 
these disturbances (see Table 3), the system was marginally sensitive. The maximum value of % 
was 0.6% which was attained after some 200 seconds had elapsed. The nonuniformity 
subsequently decayed to less than 0.2% after 1500 seconds. In order to asses the impact of the 
two thruster firing events we ran the simulation with only the series (a). While the velocity field 
responded directly to the thruster sequence, the composition nonuniformity exhibited the same 
behavior in both cases. In order to check the sensitivity limits we increased the all the amplitudes 
series (a) by a factor of ten. Again, the tenfold increase in the acceleration amplitudes resulted in a 
corresponding increase in the maximum value of i;~60%. 
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3.4. Discussion 

We have investigated the sensitivity of the Bridgman-Stockbarger technique to residual 
accelerations corresponding to single and multi-component oscillatory disturbances. The multi- 
component disturbances included simple 3-component accelerations and complex acceleration 
profiles synthesized from Fourier series analyses of acceleration data obtained on SL-3. For the 
Ge:Ga system examined, we found that for accelerations less than KL 2 Hz large amplitude 
fluctuations were found to occur in ^ for acceleration amplitudes above 10' 6 g. The system 
response to 10' 5 g accelerations at 10* 3 Hz led to composition non-uniformities on the order of 
60%. For single component disturbances, larger accelerations are required to give significant non- 
uniformities as the frequency is increased. The most noticeable difference in the response of the 
Ge:Ga system to different frequency accelerations was in the behavior of the velocity and 
composition transients as the system adjusted from the zero velocity initial state to the final state 
where the system oscillated with the frequency of the driving force. The difference is clearly seen 
upon comparison of the results for 10' 3 with those of higher frequencies. The maxium velocity 
behavior, for example, shows a distinct transient behavior for frequencies greater than or equal to 
10’ 2 Hz, which is not observed at 10' 3 Hz. For the 10‘ 3 Hz frequency acceleration, the 
nonuniformity exhibits large amplitude oscillations with twice the frequency of the acceleration, 
although there is a phase lag during the first two periods. For higher frequency disturbances, the 
larger accelerations required to obtain significant non-uniformities are sufficient to modify the 
solute field such that the modulation of % is barely noticeable at early times owing to the steady 
riase in the nonuniformity magnitude during the velocity transient . Between 200 and 300 seconds 
the composition nonuniformity reaches a maximum and clearly exhibits modulations with the 
frequency of the driving force. As time proceeds the solute field slowly adjusts to the regularly 
oscillating velocity field. By 2000-3000 seconds the solute field has relaxed into a reglarly 
oscillating state and the value of \ is diminished considerably. Furthermore, <; is modulated about a 
fixed non-zero value at twice the frequency of the driving force which indicates that composition 
field has finally adjusted from the state with high lateral nonuniformity which resulted from the 
velocity transient . 

The system was insenstitive to disturbances synthesized from SL-3 measurements with 
frequencies in the range 1 0 - 1 - 1 0 Hz, and only exhibited marginal sensitivity for amplitudes of the 
individual frequency components that were a factor of ten higher than the measured amplitudes. 
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3,7. TABLES 


3.6. Tables 

Table 1. Short time response of B, to single frequency disturbances. 


magnitude [g] 

frequency [ Hz] 

duration [s] 


lO- 3 

10-2 

500 

14.0 

10- 4 

io-i 

120 

1.5 

lO- 3 

io-i 

20 

5.5 

10-3 

1 

12 

0.26 

lO' 2 

1 

2 

0.5 




Table 2. Acceleration profile components used to obtain the results depicted in Fig. 20 



HHi 

PHASE 

(radians) 



PHASE 

(radians) 

0.667S974E-05 

0.322265 6E+01 

-0.1023456E+00 

0.4157672E-05 

0.6591797E+01 

0.1069476E+01 

0.1 168701E-04 

0.3369141E+01 

-0.2171423E+01 

0.8624586E-05 

0.6738281E+01 

0.2299790E+01 

0.4029183E-05 

0.3515625E+01 

-0.2680568E+01 

0.486001 9E-05 

0.6884766E+01 

-0.5513560E+00 

0.576761 IE-05 

0.3662109E+01 

-0.1120755E+01 

0.1024745E-04 

0.7031250E+01 

0. 144197 2E+01 

0.7336796E-05 

0.3808594E+01 

-0.2688848E+01 

0.8487802E-05 

0.7177734E+01 

-0.3315119E+00 

0.2979588E-05 

0.3955078E+01 

-0.1051489E+00 

0.9614520E-05 

0.7324219E+01 

-0.2408829E+01 

0.7084803E-05 

0.4101563E+01 

-0.1396236E+01 

0.1111080E-04 

0.7470703E+01 

-0.1942924E+01 

0.2380885E-05 

0.4248047E+01 

0.4626593E+00 

0.40233 13E-05 

0.7617188E+01 

0.2730639E+O1 

0.1256049E-04 

0.4394531E+01 

-0.1144225E+01 

0.1338748E-05 

0.7763672E+01 

-0.6898107E-01 

0.1908798E-04 

0.4541 01 6E+01 

-0.2954485E+01 

0.1156079E-04 

0.7910156E+01 

-0.1999173E+01 

0.1414271E-04 

0.4687500E+01 

0.1487755E+01 

0.2205 120E-05 

0.805 6641E+01 

0.9450552E+00 

0.1242372E-04 

0.4833984E+01 

-0.2697886E+01 

0.6399363E-05 

0.8203125E+01 

0.3106058E+01 

0.50441 14E-05 

0.4980469E+01 

-0.5252368E+00 

0.1163937E-04 

0.8349609E+01 

0.2999920E+01 

0.1122994E-04 

0.5126953E+01 

0.1805229E+01 

0.7708650E-05 

0.8496094E+01 

0.2374351E+01 

0.3820708E-05 

0.5273438E+01 

-0.5802202E+00 

0.4804912E-05 

0.8642578E+01 

0.2401585E+01 

0.6576809E-05 

0.5419922E+01 

0.1198942E+01 

0.8501 992E-05 

0.8789063E+01 

0.2369910E+01 

0.3735003E-05 

0.5566406E+01 

0.1164908E+00 

0.1167328E-04 

0.8935547E+01 

0.2706169E+01 

0.2560308E-05 

0.5712891E+01 

0.2914079E+01 

0.9668714E-05 

0.9082031E+01 

0.2292677E+01 

0.2788814E-06 

0.5859375E+01 

0.1549920E+01 

0.5459605E-05 

0.9228516E+01 

0.2345090E+01 

0.1614620E-04 

0.6005859E+01 

-0.1767729E+01 

0.9634379E-05 

0.9375000E+01 

0.1297546E+01 

0.1469397E-04 

0.6152344E+01 

-0.2948359E+01 

0.4274781E-05 

0.9521484E+01 

0.2430201E+01 

0.6047145E-05 

0.6298 828E+01 

-0.1647892E+01 

0.7036096E-05 

0.9814453E+01 

0.3016870E+01 

0.2211998E-04 

0.6445313E+01 

0.2799345E+01 

0.4917858E-05 

0.9960938E+01 

-0.3074690E+01 
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0.205 1083E-05 


0.387701 3E-0S 


0.3535946E-05 


0.2020149E-05 


0.4987608E-05 


0.2028721E-05 


0.2897678E-05 


0.3776532E-05 


0.4340377E-06 


0.1476338E-05 


0.1008826E-05 


0.4067975E-05 


0.3092643E-05 


0.2499588E-05 


FREQUENCY 

(Hz) 


0.3662109E-01 


0.732421 9E-01 


0.109 8633E+00 


0.1464844E+00 


0.1831055E+00 


0.2197266E+00 


0.2563477E+00 


0.2929688E+00 


0.3295898E+00 


0.36621 09E+00 


0.4028320E+00 


0.4394531E+00 


0.4760742E+00 


0.5126953E+00 



-0.1539647E+00 


-0.1277694E+01 


0.2503309E+01 


-0.8924237E+00 


0.9381254E-01 


-0.6826349E+00 


-0.2601777E+01 


-0.1353911E+01 


0.1620389E+01 


Table 3 b. Acceleration components for thruster firing 



FREQUENCY 

PHASE 




(g) 

(Hz) 

(radians) 


Keel ■ 


0.1 140294E-04 

0.1464844E+00 

-0.6030687E+00 

0.7271495E-04 

0.5126953E+01 

-0.6024272E+00 | 

0.9509653E-05 

0.2929688E+00 


0.5785038E-04 

0.5273438E+01 


0.3754365E-05 

0.4394531E+00 

-0.5223967E+00 

0.5360622E-04 

0.5419922E+01 

-0.2476458E+00 

0.9918378E-05 

0.5859375E+00 

-0.2804558E+01 

0.2124587E-04 

0.5566406E+01 

0.1014042E+01 

0.2464890E-05 

0.7324219E+00 

-0.3086868E+00 

0.5234673E-04 

0.5712891E+01 

0.8560771E+00 


0.8789063E+00 

0.2448875E+01 

0.7330361E-05 

0.5859375E+01 

0.7808214E-01 

0.8610295E-05 

0.1025391E+01 

-0.249275 8E+01 

0.5052985E-04 

0.6005859E+01 

0.3010761E+01 

0.3151318E-05 

0.1171875E+01 

-0.2377229E+00 

0.4433946E-04 

0.6152344E+01 

-0.4996738E+00 

0.1032598E-04 

0.1318359E+01 

-0.7560444E+00 

0.552681 6E-04 

0.6298828E+01 

-0.2421690E+01 

0.8860652E-05 

0.1464844E+01 

-0.8895736E+00 

0.5564178E-04 

0.6445313E+01 

0,309 3696E+00 

0.1 105441E-04 

0.1611328E+01 

0.9257343E-01 

0.4753912E-04 

0.6591797E+01 

-0.1798343E+01 

0.1503079E-05 

0.1757813E+01 

0.2136720E+01 

0.8398587E-04 

0.6738281E+01 

0.1045000E+01 

0.921 6345E-05 

0.1904297E+01 

-0.2327524E+01 

0.8855044E-05 

0.6884766E+01 

-0.9388244E+00 

0.6179209E-05 

0.2050781E+01 

0.2664295E+01 

0.7329257E-04 

0.7031250E+01 

0.2976783E+01 

0.1986734E-04 

0.2197266E+01 

-0.2010437E+01 

0.8345300E-04 

0.71777 34E+01 

EEBH 

0.1236275E-04 

0.2343750E+01 

0.6686478E+00 

0.8957948E-04 

0.74707 03E+01 

-0.3556173E+00 

0.3020275E-04 

0.2490234E+01 

0.2269750E+01 

0.3759252E-04 

0.76171 88E+01 

-0.1476801E+01 

0.2345829E-04 

0.2636719E+01 

0.1499383E+00 

0.4840786E-04 

0.7763672E+01 

0.1300743E+01 

0.1733246E-04 

0.2783203E+01 

-0.301 6654E+01 

0.7912557E-04 

0.7910156E+01 

-0.8301902E+00 

0.1302133E-04 

0.2929688E+01 

-0.1257968E+01 

0.7635894E-04 


0.1833709E+01 

0.6113684E-05 


0.2914706E+00 

0.4677792E-04 

0.8203125E+01 

0.904309 2E-01 

0.1626974E-04 

0.3222656E+01 

-0.2269336E+01 

0.6698716E-04 

0.8349609E+01 

-0.2321067E+01 

0.8334577E-05 

0.3369141E+01 

-0.29603 12E+01 

0.1024659E-03 

0.8496094E+01 

0.8468162E+00 

0.2110433E-04 

0.3515625E+01 

-0.1538641E+01 

0.1090984E-03 

0.8642578E+01 

-0.1963995E+01 

0.2329673E-04 

mmmgm 

0.2187939E+01 

0.7792140E-04 

0.8789063E+01 

0.1164153E+01 

0.4597104E-04 

0.3808594E+01 

-0.1761344E+01 

0.9283071E-04 

0.8935547E+01 

-0.94303 96E+00 

0.1087029E-04 

0.3955078E+01 

0.5420735E+00 

0.1097537E-03 

0.908 2031E+01 

0.2549556E+01 

0.2991441E-04 

0.4101563E+01 

-0.1586179E+01 

0.1056685E-03 

0.9228516E+01 

-0.4995980E+00 

0.1452386E-04 

0.4248047E+01 

-0.44541 26E+00 

0.1023078E-03 

0.9375000E+01 

-0.2675524E+01 

0.8852465E-05 

0.4394531E+01 

0.2433605E+01 

0.701 2900E-04 

0.9521484E+01 

0.1263102E+01 

0.3346926E-04 

0.4541016E+01 

-0.6390298E+00 

0.2513885E-03 

0.9667969E+01 

-0.1458808E+01 

0.3887603E-04 

0.4687500E+01 

0.2700060E+01 

0.2497189E-03 

0.9814453E+01 

0.14469 93E+01 

0.6244 237E-04 

0.4833984E+01 

0.901 49 22E +00 

0.3261 188E-03 

0.9960938E+01 

-0.1155957E+01 

0.300641 8E-04 

0.49804 69E+01 

-0.3005708E+01 
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3.7. Figures and Figure Captions 


Fig. 1 

Fig. 2 
Fig. 3 
Fig. 4 

Fig. 5 
Fig. 6 

Fig. 7 

Fig. 8 

Fig. 10 

Fig. 11 
Fig. 12 

Fig. 13 
Fig. 14 
Fig. 15 

Fig. 16 


Maximum u-velocity a) and v-veiocity b) as a function of time for a sinusoidal 
acceleration with magnitude 10‘ 5 g and frequency 10" 3 Hz oriented parallel to the 
crystal surface. 

Compositional nonuniformity ^ as a function of time for a sinusoidal acceleration 
with magnitude 10' 5 g and frequency 10' 3 Hz oriented parallel to the crystal 
surface 

The solute field at a) 500, b) 1000, c) 1500 and d) 2000 seconds for a sinusoidal 
acceleration with magnitude 10* 5 g and frequency 10' 3 Hz oriented parallel to the 
crystal surface. 

Maximum u-velocity a) and v-velocity b) as a function of time for a sinusoidal 
acceleration with magnitude 10‘ 4 g and frequency 10 -2 Hz oriented parallel to the 
crystal surface. 

Compositional nonuniformity \ as a function of time for a sinusoidal acceleration 
with magnitude 10' 4 g and frequency 10‘ 2 Hz oriented parallel to the crystal 
surface. 

Maximum u-velocity a) and v-velocity b) as a function of time for a sinusoidal 
acceleration with magnitude 10' 3 g and frequency 10' 1 Hz oriented parallel to the 
crystal surface. 

Compositional nonuniformity \ as a function of time for a sinusoidal acceleration 
with magnitude 10‘ 3 g and frequency 10' 1 Hz oriented parallel to the crystal 
surface. 

Velocity fields at a) 0.25 b) 0.5 c) 0.75 and d) 1 second caused by a sinusoidal 
acceleration with magnitude 10' 2 g and frequency 1 Hz oriented parallel to the 
crystal surface. 

Maximum u-velocity a) and v-velocity b) as a function of time from 802 to 808 
seconds for a sinusoidal acceleration with magnitude 10' 2 g and frequency 1 Hz 
oriented parallel to the crystal surface. 

The solute field at a) 1, b) 6 c) 200 and d) 800 seconds for a sinusoidal acceleration 
with magnitude 10' 2 g and frequency 1 Hz oriented parallel to the crystal surface. 
Compositional nonuniformity \ as a function of time for a sinusoidal acceleration 
with magnitude 10‘ 2 g and frequency 1 Hz oriented parallel to the crystal surface. 
Note the low amplitude modulation of \ with a frequency of 1 Hz between 802 to 

808. Beyond 1200 seconds % fluctuates with a frequency of 2 Hz. 

The velocity fields at a) 3 b) 6 c) 8 and d) 10 seconds caused by a sinusoidal 3- 
component acceleration with magnitudes 10' 4 g, 10" 3 g and 10' 2 g frequencies 10‘ 2 , 
10' 1 and 1 Hz oriented parallel to the crystal surface. 

Maximum u-velocity a) and v-velocity b) as a function of time for a sinusoidal 3- 
component acceleration with magnitudes 10' 4 g, 10' 3 g and 10' 2 g frequencies 10' 2 , 
10" 1 and 1 Hz oriented parallel to the crystal surface. 

The solute field at a) 3, b) 6 c) 8 and d) 10 seconds for a sinusoidal 3-component 
acceleration with magnitudes 10" 4 g, 10" 3 g and 10' 2 g frequencies 10‘ 2 , 10' 1 and 1 

Hz oriented parallel to the crystal surface. 

Maximum u-velocity a) and v-velocity b) from 0-10 seconds as a function of time 
for a sinusoidal 3-component acceleration with magnitudes 10'^g, lO^g and 10‘ 3 g 
frequencies 10“ 2 , 10' 1 and 1 Hz oriented parallel to the crystal surface. 
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Fig 17 

Fig. 18 

Fig. 19 

Fig. 20 
Fig. 21 
Fig. 22 

Fig. 23 

Fig.24 


Maximum u-velocity a) and v-velocity b) for 0-200 seconds as a function of time 
for a sinusoidal 3-component acceleration with magnitudes 10' 5 g, lO^g and 10' 3 g 
frequencies 10' 2 , 10 _1 and 1 Hz oriented parallel to the crystal surface. 
Compositional nonuniformity £ as a function of time for a sinusoidal 3-component 
acceleration with magnitudes 10" 5 g, lO^g and 10* 3 g frequencies 10" 2 , 10* 1 and 1 
Hz oriented parallel to the crystal surface, a) 0 -1200 seconds, b) 1200-5500 

seconds. Note that beyond 2000 seconds, % can be seen to fluctuate with a 
frequency of 2 Hz. 

The velocity components at a point in the melt as a function of time for a sinusoidal 
3-component acceleration with magnitudes 10' 5 g, 10' 4 g and 10' 3 g frequencies 10' 
2 , 10' 1 and 1 Hz oriented parallel to the crystal surface. 

Residual acceleration profile with the amplitudes and frequencies listed in Table 2 
taken from acceleration data obtained on Spacelab 3. 

The velocity components at a point in the melt as a function of time for the 
acceleration profile shown in Fig. 20 oriented parallel to the crystal surface. 
Maximum u-velocity a) and v-velocity b) for 0-200 seconds as a function of time 
for an acceleration with the frequencies listed in Table 2 and amplitudes multiplied 
by a factor of 10. The acceleration is oriented parallel to the crystal surface. 

Compositional nonuniformity % as a function of time for an acceleration with the 
frequencies listed in Table 2 and amplitudes multiplied by a factor of 10. The 
acceleration is oriented parallel to the crystal surface. 

Lateral non-uniformity £ as a function of time for accelerations taken from a sample 
time series (see inset) constructed from data obtained on Spacelab 3. The 
acceleration consists of a repeated “noise” segment (length 10 seconds, 10' 1 Hz </ 
< 1 Hz) and a thruster firing event (length 10 seconds, 10' 1 Hz </< 13 Hz). The 
latter is introduced at 10 and 80 seconds. 
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FIG. 13 
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4. 


MODELLING CRYSTAL GROWTH BY PHYSICAL VAPOR 
TRANSPORT IN CLOSED AMPOULES 


4.1. Summary 

There has been great interest in growing crystals by physical vapor transport in closed 
ampoules ever since mercuric iodide and cadmium telluride crystals were successfully grown on 
Spacelab missions in 1985. In order to ascertain if the transport modes prevailing in low gravity 
environments are beneficial to vapor growth and to determine tolerance limits of residual 
accelerations for such growth, we have developed a steady state two-dimensional model. This was 
solved using the finite volume code PHOENICS. Gravitational accelerations of less than O.lg 
were found to be sufficient to suppress bouyancy-driven convection to an extent that diffusion was 
the dominant transport mode. Furthermore, we found that the convection is usually solutally- 
driven and that if convection is reduced enough there is greater uniformity in the growth rate. It is 
also shown that a linear temperature gradient on the ampoule walls will cause the vapor to be 
supersaturated throughout the ampoule resulting in undesireable nucleation at the walls. A 
"hump" in the wall temperature profile can be used to avoid this. The size of the hump could be 
used to avoid this. The prevailing transport conditions determine the size of the hump needed. 


4.2. Introduction 

In 1985 van den Berg and co-workers grew a single crystal of mercuric iodide (Hgl 2 ) on 
Spacelab 3 by physical vapor transport [1]. That crystal was shown to have significantly better 
structural and electronic properties than Earth-grown crystals. Two reasons were suggested for this 
improvement in quality. The first one is the propensity for slippage of Hgl 2 crystals under 
terrestrial gravity and the resulting decrease in structural quality. The second is the reduction of 
convective flows in a low gravity environment which was expected to produce a more uniform 
growth rate and, hence, crystals of greater homogeneity. Even if the fonner reason was specific to 
soft crystals like HgU, he latter, if correct, would suggest that a low gravity environment is 
beneficial for all crystals grown by physical vapor transport (PVT). These improvements in crystal 
quality realtive to earth-grown crystals were further confirmed by Nitsche and co-workers [2] for 
the vapor growth of cadmium telluride crystals on the Spacelab D-l mission in 1985. 

These results and those of the earlier experiments with GeSe and GeTe by Wiedemeier and 
co-workers [3,4] have sparked interest in growing other crystals by PVT in future shuttle 
missions and on the Space Station. The interplay of diffusion and convection in such systems is 
very complex and quantitative results cannot be predicted by simple models or purely physical 
arguments. However, as has been demostrated for crystal growth from the melt and solution, 
"generic" models of the growth process can be used to identify trends in the response of the 
condition to different transport conditions. In this work we will attempt to arrive at a better 
understanding of the role of convection in a simplified PVT model and to predict the effect of low 
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gravity on the growth system. In addition, it is hoped that the results of this study will lead to 
improvements in the growth methods employed. 

The first systematic attempt to model PVT in a closed ampoule with convective flow was 
made by Klosse and Ullersma [5]. They ignored the net mass transport resulting from crystal 
growth and also solutal convection which makes their results of limited validity in real growth 
situations. Rosenberger and co-workers have demonstarted the complexity of the transport even in 
a simplified two-dimensional model, showing that recirculation of the inert gas can result in 
convective flows even in the absence of gravity [6]. Also they incorporated the compressible form 
of the Navier-Stokes equation (rather than the Boussinesq approximation), used the exact 
boundary conditions at the source and crystal interfaces to obtain the most realistic model used so 
far for PVT [7-9]. They were able to investigate the effect of a variety of growth parameters like 
aspect ratio and the orientation of the ampoule with respect to gravity. More importantly, they 
made clear the nature of the flow field during PVT and showed that non-uniformities in the growth 
rate can occur even in the best of situations. 

Bontoux and co-workers have modelled 3-D convective flow during PVT [10,11]. Despite 
their use of the Boussinesq approximation their results compared favourably with experiments 
[12]. The relationship of the crystal growth rate to the convective flow for different experimental 
conditions cannot be readily deduced from their work. 

In our study we will concentrate mainly on exploring the advantages or disadvantages of 
low gravity environments for vapor growth though some of the conclusions should also be 
relevant to ground-based experiments. We will employ a "generic" model which includes aspects 
of PVT important to crystal growers which have been neglected in previous models, such as the 
use of nonlinear wall temperature gradients to prevent spurious nucleation, and the use of different 
total pressures and inert gases in the ampoule. We will present our model in the next section and 
the numerical scheme employed in section 3. In section 4 present and discuss our results and draw 
conclusions in section 5. 

4.3. The Model 

We restrict our model to examine how a typical convective flow field will affect the crystal 
growth rate. A two-dimensional model will suffice and we have used a rectangle to represent the 
ampoule with a source and crystal of material A on the two ends. In addition, the domain contains 
an inert gas B. 

The other major assumption we have made is to postulate the existence of a basic steady 
state. While this can never be achieved in real growth situations, it is nonetheless usfeul for a 
qualitative interpretation of the response of the system to different transport conditions. Not only is 
the temperature gradients at the walls hard to maintain with moving source and crystal interfaces 
but, more importantly, a steady state would require dopant incorporation at the same rate as 
sublimation at all times, a very unlikely prospect. Time-dependent calculations are a time 
consuming and unless absolutely necessary should be avoided. Besides, the main interest in using 
a time dependent model is to study dopant incorporation and compositional uniformity, but this 
also requires a knowledge of the composition of the vapor and the vapor-solid equilibrium 
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relationships for the dopant and host crystal. Unfortunately for most systems this information is 
not available. Thus, the effort involved in such detailed calculations is of limited practical use. 

We have also assumed that all physical properties except density and diffusivity are 
constant and this is reasonable for most vapor growth systems. The conservation equations for a 
compressible gas then become: 



Vpv = 0, 

(1) 


pv-Vv = -Vp + pV 2 v + pg, 

(2) 


v-VT = kV 2 T, 

(3) 

and 

v-VW a = V D AB VW a . 

(4) 


Here v, p, T, Wa and g are the mass averaged velocity [13], pressure, temperature, mass 
fraction of the transported material A and gravitational acceleration respectively. The physical 
properties p, (X, k and Dab are the density, viscosity, thermal diffusivity and the binary diffusion 
coefficient, respectivley. We have assumed an ideal gas relationship in relating the density to the 
pressure. 

The boundary conditions on velocity are that the mass averaged velocity vanishes at the 
ampoule walls (for a discussion of this point see [6,14]) and that at the source and crystal ends 
there will be a perpendicular component of velocity due to the inflow and outflow of component A. 
These conditions are 

v = 0 at x = 0 and x = a, (5) 


vn = - VWa-ii at y = 0 and y = l. 

1-Wa (6) 

The two ends are held at different constant temperatures. The walls are maintained at a 

specified temperature profile. We will not consider the case of insulated walls as this is difficult to 

achieve in practice and also since most vapor growth experiments are conducted with an imposed 


temperature profile to avoid nucleation at the ampoule walls. Thus, 

T = T s at y = 0 and T = T c at y = l, (7) 

and T = f(y) at x = 0 and x = a . (8) 

The boundary conditions on the mass fraction of component A are impermeability at the 
ampoule walls and vapor-solid equilibrium at the source and crystal: 

VWaii =0 at x = 0 and x = a ( 9 ) 

and p A = Peq(T) at y = 0 and y = /, (10) 


where p A is the partial pressure of component A and Peq(T) is the equilibrium vapor pressure of A. 
Since we will be using the iodine-hydrogen system for most of our calculations, with iodine the 
crystallizing material and hydrogen the inert gas, we have used the relationship suggested by 
Gillespie and Fraser [15] for the vapor pressure of iodine: 

P eq (T) = 760 exp( 30.795 - 8088.64/T - 2.0131nT ) mm Hg. (11) 

The diffusion coefficient is a strong function of temperature and pressure but for vapor 
growth the total pressure is essentially constant over the ampoule and the temperature differences 
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are usually too small to cause significant variations. Nevertheless some simulations were carried 
out where the diffusion coefficient, as suggested by Rosenberger [13], was given by 

Dab = D 0 (p 0 /p)(TA’o) 1,82 » (12) 

where Do is the diffusion coefficient at temperature To and pressure p 0 . However, the results so 
obtained hardly differed from those with constant coefficients. Hence, most of the simulations 
were carried out with constant coefficients though for each simulation at an assumed total pressure 
Dab was recalculated using eqn. (12). The physical properties used in the simulations are listed in 
table 1. 


4.4. Numerical Method 

The coupled conservation equations (1-4) were solved subject to the boundary conditions 
(5-11) using the commercial code PHOENICS. This code is based on a finite volume method with 
a fully implicit time scheme and allows the solution of the transport equations in the compressible 
form [16]. Convergence was not easy to obtain and the mesh size and the relaxation parameters 
had to be adjusted for each simulation. A particular difficulty was to have mass conservation 
between inlet and outlet, as this is not automatically satisfied when the boundary conditions are 
discretized. The mesh we used was finer near the crystal and the source and coarser elsewhere. 
Typically a simulation with an ampoule of aspect ratio (Length/Width) //a = 2 required a 24x36 
mesh and one with an aspect ratio of 5 required a 24x52 mesh. 

All computations were carried out on the Cray-XMP of the Alabama Supercomputer 
Network. Computer time for each run varied according to the complexity of the formulation (eg. 
the use of a polynomial for the imposed wall temperature profile and the use of a temperature 
dependent diffusion coefficient), the magnitude and orientation of the gravitational acceleration and 
the relative magnitudes of the thermal and solutal density gradients when they opposed each other. 
Most of the simulations required between 5 to 10 minutes of Cray CPU time. 

4.5. Results and Discussion 

The first question examined was the effect of gravitational acceleration on the crystal 
growth rate. Greenwell et al. [8] first demonstrated that the convective flow can cause 
nonuniformities in the growth rate. We examined this in more detail using a linear temperature 
profile at the ampoule walls. It was found that when solutal convection dominates, the imposed 
temperature profile has no effect on the growth rate. For most of our calculations weused the 
properties of iodine as the crystallizing material A and hydrogen as the inert gas B. For some of 
the calculations we varied the molecular weight of B from 2 up to 254 in order to study the effect 
of thermally driven convection in isolation. 

Fig. 2 shows our results for a horizontal ampoule of aspect ratio 5 (length=5cm, 
width=lcm)), with a temperature of 370.5 K at the source end and 358.1 K at the crystal end, and 
a total pressure of 100 torr. The most important observation from the point of view of space 

processing is that even at 0.1 g convective effects become very small and we have essentially 
diffusive transport in the vapor. Moreover, when the same ampoule was kept vertical (with source 
both up and down) even a gravitational level of lg was not enough to cause significant changes in 
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growth rate compared to the diffusive case. This is quite unlike other crystal growth methods 
being tested in low gravity environments which require gravitational levels of 10 _6 g or lower to 
avoid convective effects [17,18], which may be difficult to obtain in the currently available low 
gravity environments. For vapor growth these environments would seem quite adequate and 
problems of "g-jitter" which are a cause of concern for other experiments should have little effect 
on vapor growth systems. The experiments of Launay and co-workers [19] who grew Ge by 
chemical vapor transport on the Spacelab D1 mission in 1985 confirms this, as their measurements 
of the mass transport in the ampoule indicated a purely diffusive regime. However it must be 
mentioned that in their numerical model of this experiment they reported significant convection 
even at 10' 2 g and observed a purely diffusive solution only at 10 4 g [20], 

This leaves us with the question of whether there is anything to be gained by suppressing 
convection. From fig. 2 we can see that convection does enhance the mass transport of component 
A, and consequently also the growth rate, but at the expense of nonuniformity in the growth rate. 
It also can be seen that convection and, hence, the nonuniformity, is more likely to be significant 
when the ampoule is horizontal rather than vertical. In the vertical position the thermal and solutal 
gradients oppose each other and for convection to occur the solutal gradient must exceed the 
thermal gradient and overcome the viscous resistence. In the horizontal position convective effects 
are a more likely occurrence. In the vertical position we can have diffusive growth, especially if 
the ampoule is held in the convectively stable position. Suppressing convection should be possible 
even without resorting to a low gravity environment. 

Fig. 3 shows the results for an aspect ratio 2 (length=2cm). The growth rates are much 
larger than for aspect ratio 5 but the nonuniformity due to convection is also far greater. The 
increase in the convection level is visible from the fact that even in the vertical position convection 
will be significant. However, even for this case, a factor of ten reduction in the gravitational level 
is enough to suppress convective effects on growth rate. Increasing the temperature difference 
between source and crystal gave similar results (fig. 4). This suggests that one benefit of the low 
gravity environment would be to increase the growth rates by using ampoules of smaller aspect 
ratios or using larger temperature differences between source and crystal and still avoid the 
drawback of increased convection, which normally would rule out such methods on Earth. 

Next we examined situations for which thermal convection might be larger than the solutal. 
This could occur when inert gases other than hydrogen and having a molecular weight closer to 
that of the crystallizing material are used. Thermal convection, if dominant, is very hard to 
eliminate on Earth. Unlike the mass fraction it is difficult to maintain a uniform temperature across 
a cross section of the ampoule, which means radial temperature gradients are invariably present. 
The results for a horizontal ampoule are shown in fig. 5. As the molecular' weight of B is increased 
we can see that the effect of solutal convection first increases and then decreases rapidly and 
eventually the transport becomes largely diffusive. Its only when the molecular weights of A and 
B are almost equal that the thermal gradient overcomes the stabilizing solutal gradient. When the 
molecular weights were equal the thermal gradient gave rise two convection rolls, unlike the single 
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rolls of solutally driven convection seen so far, and resulted in significant nonuniformities in the 
growth rate. From this it is possible to say that unless the two components have molecular weights 
very close to each other, the dominant mode of convection is likely to be solutal. But it should be 
mentioned that when the two gradients were comparable it was difficult to get the computer code to 
converge possibly suggesting that oscillatory convection may have been present. 

A very important concern in vapor growth is the avoidance of parasitic nucleation at the 
ampoule walls [21]. Figs. 6 and 7 show the equilibrium vapor pressure of component A computed 
from eqn. (11) and the partial pressures of A at the ampoule walls resulting from the transport, 
when a linear temperature profile is used. The results are for a source temperature of 378. 1 K and 
the crystal at 358.1 K with fig. 6 for the case of purely diffusive transport and fig. 7 for the 
ampoule kept horizontal with respect to gravity. It is easy to see why a linear temperature profile is 
rarely used in practice, as the vapor of component A is in a supersaturated state throughout the 
ampoule. 

The commonly used experimental technique to prevent nucleation at the walls is to impose a 
wall temperature profile with a maximum between the crystal and the source, usually referred to as 
a temperature "hump" [21-23]. We ran some simulations in a horizontal ampoule for a source 
temperature of 378.1 K and the crystal at 358.1 K with a temperature hump of 4 K, the functional 
form of the temperature profile being: 

f(b) = 378.1 + 26.67b - 55.66b 2 + 28.04b 3 + 115.46b 4 - 134.61b 5 K (13) 

where b is y //. Fig. 8 shows the results and we can see that even though the supersaturation has 
been eliminated, the profile is an overcorrection for the problem of nucleation at the walls. A 
smaller hump would have sufficed. But this result is valid only when the diffusion coefficient is 
large as in hydrogen pairs. When inert gases of larger molecular weight are used Dab is likely to 
be much smaller. Even for hydrogen pairs at pressures close to 1 atm. Dab will be 0.16 from eqn. 
(12). In fig. 9 we have plotted the results with a diffusion coefficient of 0.25 and it can be seen 
that component A becomes supersaturated near the crystal and a larger hump would be needed to 
avoid it. This also explains why nucleation is more likely to occur near the crystal end and also 
why the temperature maximum has to be located closer to this end. It is worth mentioning that this 
problem at low diffusion coefficients occurs only due to convection. At low gravity smaller humps 
would be quite adequate to avoid supersaturation. If the variation of the diffusion coefficient with 
temperature is significant it is likely the effect will be felt most in determining the magnitude of the 
hump. We performed one simulation with Dab as a function of temperature given by eqn. (12). 
Owing to the small temperature differences employed, the results were indistinguishable from 
those for a constant coefficient . 

One consequence of using temperature humps is that when thermal convection is dominant 
it is confined to a small region near the crystal end of the ampoule as shown in fig. 10. Here the 
convection is strong enough to cause nonuniformities in the growth rate but as it is confined to a 
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small region we do not have the benefit of enhanced mass transfer across the ampoule for, as can 
be seen in fig. 1 1, there is no overall increase in the growth rate. Barring this instance, when the 
molecular weights of both components are comparable, the temperature hump had no effect on the 
growth rate. So any conclusions obtained from a linear temperature profile would be applicable 
here as well. 

4.6. Conclusions 

Convective effects can easily be suppressed in vapor growth systems for most low gravity 
environments. A gravitational acceleration level of less than O.lg is adequate to ensure purely 
diffusive transport. The more important question is then: What are the benefits of purely diffusive 
transport conditions? 

The benefits are by no means overwhelming but they do exist. Nonuniformities in the 
growth rate caused by convection can be eliminated, freeing the experimenter to use methods to 
enhance the growth rate like smaller aspect ratios and larger temperature differences between 
source and crystal, which would not have been feasible otherwise. In addition in order to avoid 
nucleation at the walls the size of the temperature hump for purely diffusive conditions can be 
smaller than that needed in the presence of than convection. 

Finally we must reemphasize that a thorough study of physical vapor transport in closed 
ampoules would require (1) an experimental program to determine the composition of the vapor 
and equilibrium relationships between vapor and solid for both host and dopant components, and 

(2) a time-dependent numerical model to investigate the effects of experimental conditions on the 
compositional uniformity of the crystals grown. Proper evaluation of the (improved) quality of 
crystals for device applications can only be made if there is a well-defined link between transport 
conditions and the characteristics of the grown crystal which are necessary for high quality device 
fabrication. 
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4.8. Tables 


Table 1. Physical properties and operating conditions 


Property 

Value 

Viscosity (fJ.) 

: ■ ... ~ ' T/. ■' ; ' 

Diffusion Coefficient 

1.25 cm^ sec -1 

at 360 K and 100 torr 


Prandlt number ((i/pK) 

1.0 

Molecular weight of A (Ma) 

254 

Operating Conditions 

Total Pressure (p) 

100 mm Hg 

Crystal Temperature (T c ) 

358.1 K 

Source Temperature (T s ) 

370.5 and 378.1 K 

Ampoule Width (a) 

1 cm 

Ampoule Length (/ ) 

2 and 5 cm 
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4.9. Figure Captions 

Fig. 1 Ampoule geometry used in the model showing location of source and crystal. 

Fig. 2 Crystal growth rates for an ampoule of aspect ratio 5 (length = 5cm) and a linear 
temperature gradient of 2.48 K cnr 1 along the walls. 

Fig.3 Crystal growth rates for an aspect ratio 2 (length = 2cm) and a linear temperature gradient 
of 6.2 K cm' 1 along the walls. 

Fig. 4 Crystal growth rates for an aspect ratio 5 and a linear temperature gradient of 4 K cm' 1 
along the walls. 

Fig. 5 Crystal growth rates for various species of B with differing molecular weights for a 
horizontal ampoule of aspect ratio 5 and a linear temperature gradient of 4 K cm' 1 along 
the walls. 

Fig. 6 Plot of the equilibrium vapour pressure and the component A partial pressure at the walls 
for diffusive transport for an ampoule with an aspect ratio of 5 and a linear temperature 
gradient of 4 K cm' 1 along the walls. 

Fig. 7 Plot of the equilibrium vapour pressure and the partial pressures of A at the walls 
horizontal ampoule of aspect ratio 5 and a linear temperature gradient of 4 K cm' 1 along 
the walls. 

Fig. 8 Plot of the equilibrium vapour pressure and the partial pressures of A at the walls when 
the ampoule is kept horizontal. The aspect ratio is 5 and the imposed wall temperature 
profile is given by (13). The temperature hump is 4K. 

.Fig. 9 Plot of the equilibrium vapour pressure and the partial pressures of A at the walls when 
the ampoule is kept horizontal for a Dab of 0.25. The aspect ratio is 5 and the imposed 
wall temperature profile is non-linear. 

Fig. 10 Velocity in a horizontal ampoule of aspect ratio 5 with a non-linear imposed wall 
temeperature profile. The molecular weights of A and B are equal. Only half the ampoule 
is shown. 

Fig. 1 1 Crystal growth rates with a non-linear temperature profile at the walls and when the 
molecular weights of components A and B are equal. The aspect ratio is 5. 
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The effects of steady and impulse-type residual accelerations on dopant distributions during directional solidification in 2D and 
3D “generic” models of the Bridgman- Stockbarger technique have been investigated using numerical models. The calculations are 
based on the thermo-physical properties of molten germanium doped with a low concentration of gallium. A novel Chebyshev 
collocation pseudo-spectral method was used for the solution of the governing momentum-, mass-, species- and heat-transfer 
equations. Only convection caused by temperature gradients was considered. It was found that lateral non-uniformity in composition 
is very sensitive to the orientation of the steady component of the residual gravity vector and to the particular operating conditions 
under consideration. For growth rates of several microns per second a steady background level of 10 1 * * * * 6 times that of normal gravity 
can be tolerated provided that the acceleration vector is aligned with the axis of the growth ampoule. For reduced growth rates, 
higher (steady) background acceleration levels are acceptable. It was also found that laterally or radially averaged composition 
profiles are alone insufficient to describe the extent of residual convection in a spacecraft environment. The effects of impulse-type 
disturbances can be severe and can extend for times on the order of one thousand seconds after the termination of the impulse. A 
so-called “compensating double pulse” will not result in completely offsetting effects. 


1. Introduction 

During the past decade there has been consid- 
erable interest in the possibility of using the al- 
most weightless environment of an orbiting 
spacecraft as a laboratory in which to study crystal 
growth processes. The attraction of such a labora- 

tory is that buoyancy-driven convection in nutri- 
ent phases may be reduced or practically elimi- 
nated so that mass transfer is more likely to be 

controlled by diffusion [1-3], It has been re- 

cognized for some time that convection can lead 
to undesirable compositional variations in melt 
grown (semi-conductor) crystals [2,4-6], In many 
of these cases the convection is buoyancy-driven 
[7]. Thus, as far as crystal growth is concerned, the 
diffusion-controlled conditions which may result 
in a spacecraft should favor compositional uni- 
formity in the crystal [1,2], 

That the magnitude of the effective gravity 
vector may be reduced as much as six orders of 
magnitude is well known, and has led to the use of 


the term microgravity to describe the spacecraft’s 
residual acceleration environment. However, our 
experience in space laboratories has led to the 
realization that, while in principle the accelera- 
tions could be as low as 10 ~ 5 -10" 6 gal *, in 
practice this is not always the case [8,9]. The many 
sources of residual acceleration [8-10] include the 
earth’s gravity gradient, atmospheric drag, basic 
orbital attitude motions, crew motions, machinery 
induced perturbations, and attitude corrections. 
The latter three can give rise to disturbances on 
the order of 10“ 3 to 10" 2 gal. While these acceler- 
ations may be orders of magnitude lower than that 
experienced at the earth’s surface, they are non- 
etheless finite and pose potential problems for 
certain types of experiment, especially those for 
which minimization of accelerations is desirable. 

In this work, we examine the extent to which 
buoyancy-driven convection, caused by small 

* 1 gal = 9.8 m s -2 , and is a unit of acceleration (named after 

Galileo) commonly used in geophysics. 
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(North-Holland Physics Publishing Division) 



286 


J.I.D. Alexander et al. / Low gravity tolerance of Bridgman-Stockbarger growth. I 


accelerations characteristic of the steady compo- 
nent of the residual acceleration vector in a 
spacecraft, result in non-uniform solute distribu- 
tion in a directionally solidified crystal. The physi- 
cal system under consideration involves crystals 
grown from a two component melt via the Bridg- 
man-Stockbarger technique. In particular, the 
importance of ampoule orientation relative to the 
gravity vector is studied. The effects of impulses 
are also examined. The consequences of other 
types of time-dependent disturbances will be dis- 
cussed in a companion paper. 

In previous work involving the use of numerical 
models to study melt convection, the thermal 
boundary conditions have been modelled with 
varying degrees of fidelity. Early work [11-13] 
includes a variety of imposed temperature bound- 
ary conditions. These range from purely vertical 
temperature gradients in which the flow ensues 
after a critical value of the Rayleigh number is 
exceeded [14], to idealized conditions associated 
with Bridgman- Stockbarger furnaces [13] which 
are imposed directly on the melt and crystal 
without consideration of the heat transfer between 
the ampoule, furnace and sample. For these 
boundary conditions, flow always occurs owing to 
the presence of radial temperature gradients. Later 
models [15,16] have accounted for the presence of 
the ampoule, and the furnace design. In an actual 
growth situation the characteristic thermal profile 
of the inner surface of the furnace is not realized 
at the ampoule wall; it is modified by heat trans- 
fer between the crystal, melt, ampoule and the 
furnace itself. The tendency is to reduce axial 
temperature gradients, while radial temperature 
gradients may increase or decrease depending on 
the specific nature of the heat transfer between the 
charge and ampoule [15]. The influence of melt 
convection on dopant distribution has been ex- 
amined for dilute and non-dilute melts [13,15]. 
For a given furnace-ampoule combination the 
amount of compositional non-uniformity (or ra- 
dial segregation) was shown to be a non-linear 
function of the Rayleigh number. The growth rate 
and physical properties of the melt also influence 
the degree of compositional uniformity in the 
grown crystal. 

That residual accelerations on the order of 


10 -6 -10~ 5 gal are sufficient to cause observable 
fluid motion has been established [10,17-19] and 
can be inferred from work such as [13]. The sensi- 
tivity to low gravity has also been estimated for a 
wide variety of experiments by dimensional analy- 
sis [20,21]. The extent to which gravity causes 
solute redistribution, via buoyancy-driven fluid 
motion, during directional solidification in geome- 
tries characteristic of the Bridgman-Stockbarger 
technique has been examined using numerical 
models [11-14,19,22], Apart from McFadden and 
Coriell [22] who examined a two-dimensional 
situation in which the gravity vector continuously 
rotates relative to the container, and Polezhaev et 
al. [19], other work is restricted to axisymmetric 
situations in which a steady gravity vector is ori- 
ented parallel to the axis of a cylindrical ampoule. 

A general conclusion that can be drawn from 
all previous attempts to characterize gravity-driven 
convective effects on directional solidification from 
two component melts (and from the results of the 
work presented here) is that the maximum lateral 
solute non-uniformity (radial segregation for the 
axisymmetric cases) occurs near the transition from 
diffusion dominated to convection dominated 
growth conditions [20,23]; that is, when convective 
velocities are of the same order of magnitude as 
the diffusive velocities. The conditions under which 
this “transition” takes place will depend on the 
specific nature of the forces driving convection. 
We demonstrate in this work that the orientation 
of the steady component of the gravity vector is 
crucial in determining the magnitude of the grav- 
ity vector at which this transition occurs. Thus, for 
a given set of operating conditions the orientation 
of the gravity vector can determine the suitability 
of a low gravity environment for directional solidi- 
fication experiments. 

In the following sections we formulate the basic 
model, define the range of operating conditions 
under consideration and examine the effects of 
magnitude and orientation of steady and impulse- 
type acceleration vectors on compositional uni- 
formity in directionally solidifying crystals. The 
objective of this work is to identify trends. Rather 
than attempt to model the thermal conditions 
corresponding to a specific furnace, we adopt a 
generic model which is based on the pioneering 
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work of Chang and Brown [13], in which the 
thermal profile is imposed on the melt and crystal. 


2. Formulation 

The following physical situation, which corre- 
sponds to directional solidification by the Bridg- 
man-Stockbarger technique, is considered. An 
ampoule containing a (dilute) two-component melt 
is translated at a constant rate between fixed hot 
and cold zones of a furnace. The upper and lower 
parts of the furnace are maintained at a uniform 
temperature by means of heat pipes and the two 
zones are separated by a thermal barrier [24], 
Directional solidification takes place as the 
ampoule is translated. In an actual growth situa- 
tion the heat transfer conditions between the 
crystal, melt and ampoule can result in a non- 
planar crystal-melt interface (for a summary, see 
ref. [24]). Curvature of the solid -liquid interface 
can result in significant lateral compositional 
non-uniformity [25,26]. Since we wish to focus 
attention on the influence of convection on the 
composition of the crystal we choose to constrain 
the interface to be planar. 

Two- and three-dimensional models are consid- 
ered. For the 2D model the dilute binary melt is 
assumed to occupy a rectangular region Q of 
height L and width W , which is bounded by 
planar surfaces (see fig. 1). In the 3D model the 
rectangular region is replaced by a circular cyl- 
inder with radius W/2. 

Translation of the ampoule is modelled by 
supplying a doped melt of dilute bulk composition 
at a constant velocity V M at the top of the 
computational space, and withdrawing a solid of 
composition c s (which, in general, will be a func- 
tion of both space and time) from the bottom. The 
crystal-melt interface is located at a fixed dis- 
tance L from the top of the computational space. 
Thus, the assumption that the ampoule translation 
rate and the growth rate are equal is implicit. The 
thermal barrier which separates the zones is mod- 
elled using adiabatic sidewalls of length L / 4 (see 
fig. 1). The temperature at the interface is taken to 
be r M , the melting temperature of the crystal, 
while the upper boundary is held at a higher 


W=2R 



Fig. 1 . Idealized model of the Bridgman-Stockbarger system. 


temperature T H . In an actual experiment, owing to 
the finite length of the ampoule there is a gradual 
decrease in length of the melt zone during growth. 
In this model transient effects related to this 
change are ignored. Thus, it is assumed that the 
ampoule is sufficiently long for steady state tem- 
perature. concentration and growth rate to be 
established. The appropriate boundary condition 
to adopt in this case is the pseudo-steady state 
condition used by Chang and Brown [13], The 
only transient effects to be considered will arise 
directly from the time-dependent nature of the 
residual gravity field. We also assume that the 
contribution of the solute (dopant) to convection 
is negligible. Free convection is driven by thermal 
gradients only. 

The governing equations are cast in dimension- 
less form using L, k/L (where k is the melt’s 
thermal diffusivity), p M ic 2 /L 2 (where p M is the 
melt’s density), T n — T m , and c x to scale the 
lengths, velocity, pressure, temperature and solute 
concentration, respectively. The dimensionless 
equations governing momentum, heat and solute 
transfer in the melt are then 

jf- + (grad u) u 

= —grad p + Pr Au + Ra Pr 8 g(t), 
div u — 0, 


( 1 ) 

( 2 ) 
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Table 1 

Forms of the acceleration vector examined in this work 
Steady 

So = So. A + go.-*! 

Ilgll = itr 4 , y^xitr 5 , ltr 5 , sfi xicr 6 , sxitr 6 , v^xio -6 , kt 6 , icr 7 

Impulse 

g(0 = [go + £(')]*> g 0 = HT 6 ,^ xl(T 6 , 

g(t) = 0 for t < t r and t > t 2 and g(t ) — gj for t\<t < t 2 

g(t) = 0 for t < q, t 2 < I < ?3 and t > t 4 and g(t ) = g l for t x <t <t 2 , t 3 <t <t 4 

g(t ) = 0 for t < / 1 , t 3 <t < t 2 and t > t 4 and g(t) = g , for q < / < t 2 , g (f) = — g, for q < t < t 4 , gj = 5xlO~ 3 , 3xl0 -3 


00 

+ u ■ grad 8 = Ad, 

(3) 

Sc / 9C \ 

(4) 

-(^ + «.gradc)=dC, 


where u(x, t), 6 = (T(x, t) - T M )/(T H - T M ) and 
C, respectively represent the velocity, temperature 
and solute concentration. The parameters Pr = 
v/k, Ra = gft(T H — T m )L 2 /kv and Sc = v/D are, 
respectively, the Prandtl, Rayleigh and Schmidt 
numbers. The term g(t ) in eq. (1) represents the 
steady or time-dependent gravity vector. The value 
of g in Ra is taken to be 1 gal, i.e. equal to the 
gravitational acceleration at the earth’s surface. 
Thus, the magnitude of g represents the ratio 
between the actual residual acceleration and g. 
Table 1 lists the forms of g(t) used in these 
calculations. 

The following boundary conditions are applied 


at the crystal-melt interface 


0 = 0, 

(5) 

u ■ N — P e/a, 

(6) 

« X TV = 0. 

(7) 

3 C Pe Sc . . 

3z P, (1 k) C ’ 

(8) 


range in concentration given by 


£ = 


(C sma x - C smin ) X 100% 


where c s is the (dimensional) solute concentration 
in the crystal, and c av is the average interface 
concentration. At the “inlet” (z = 0) the following 
boundary conditions are applied 


0 = 1 , 


« ■ N = Pe a , 
« X 7V= 0. 


9 C 
9 z 


Pe Sc 
Pr 


(c-i), 


(9) 

( 10 ) 

( 11 ) 

( 12 ) 


Eqs. (8) and (12) express conservation of mass at 
the crystal-melt interface and the inlet respec- 
tively. Eqs. (6) and (10) guarantee continuity of 
the melt with the crystal and with the supply of 
melt at the inlet, while eqs. (7) and (11) ensure 
no-slip tangent to the interface and the top surface. 
At the side walls the following conditions are 
applied 

gradC-e w = 0, i/-iV=Pe o, e w • u = 0, (13) 

along with 


where N points into the melt and is the unit 
vector perpendicular to the planar crystal melt 
interface, Pe = V ml /k is the Peclet number, a = 
Pm/Ps and k is the segregation coefficient. We 
define the measure of compositional non-uniform- 
ity in the crystal at the interface to be the lateral 


0=1 (14) 

in the isothermal zone, and 
grad 0 • e w = 0 (15) 

in the adiabatic zone. Here e w is the normal to the 
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Table 2 

Thermo-physical properties of gallium-doped germanium [13] 
and operating conditions used in the calculations 


Property 

Value 

Thermal conductivity of the melt 

0.17 W K -1 cm -1 

Heat capacity of the melt 

0.39 Jg” 1 K" 1 

Density of the melt (p M ) 

5.6 g cm -3 

Density of the solid (p s ) 

5.6 g cm -3 

Kinematic viscosity of the melt (r) 

1.3XlO~ 3 cm 2 s -1 

Melting temperature ( T M ) 

1231 K 

Solute diffusivity ( D ) 

1.3xl0“ 4 cm 2 s ~ 1 

Thermal diffusivity of the melt (k) 

1.3X10~ ] cm 2 s _1 

Segregation coefficient (k) 

0.1 

Thermal expansion coefficient (/3) 

2.5X10 -4 K _1 

Associated dimensionless parameters 

Value 

Prandtl number Pr = v/k 

0.01 

Peclet number Pe = V M L/ k 


Schmidt number Sc = v /D 

10 

Density ratio a 

1.0 

Operating conditions 

Value 

(A) Hot zone temperature (T H ) 

1331 K 

Distance between inlet 


and interface (L) 

1.0 cm 

Height of adiabatic zone 

2.5 mm 

Ampoule width (diameter) 

1.0 cm 

(B) Hot zone temperature (r H ) 

1346 K 

(C) Hot zone temperature (T H ) 

1251 K 

(D) Ampoule width 

2.0 cm 

(E) Ampoule width 

0.5 cm 

(VI) Translation (supply) rate (F M ) 

6.5 jam s~ 1 

(V2) Translation (supply) rate (F M ) 

3.25 ji m s _1 

(V3) Translation (supply) rate (F M ) 

0.65 pm s _1 


ampoule wall. While the above model will rarely 
apply to a specific furnace (since details of the 
heat transfer at the ampoule walls, and between 
ampoule and furnace are neglected), it nonetheless 
serves as a reasonable “generic” model with which 
to carry out a preliminary analysis of a directional 
solidification experiment under conditions char- 
acteristic of the low gravity environment of space. 

For the most part, our calculations are limited 
to thermo-physical properties corresponding to di- 
lute gallium-doped germanium, although the effect 
of different Schmidt numbers are considered for a 
limited range of steady low gravity conditions. 
The values of the thermo-physical properties and 
the associated dimensionless groups and operating 
are given in table 2. 


3. Method of solution 

The governing equations (1)-(15) were solved 
by two different methods; a pseudo-spectral 
Chebyshev collocation method and a finite volume 
method implemented by the code PHOENICS 
[28]. The latter method was used to obtain most of 
the results pertaining to the Ge : Ga system. Only 
those results corresponding to B (2D and 3D) in 
table 2 and to the 2D impulse accelerations were 
computed using PHOENICS. 

The pseudo-spectral method employed is an 
orthogonal collocation pseudo-spectral method 
which uses Chebyshev polynomial expansions for 
spatial discretization [27,29]. Only the two-dimen- 
sional calculations were tackled using this ap- 
proach, although the method is readily extended 
to three dimensions. Steady solutions are obtained 
using a pseudo-unsteady method which involves a 
generalized ADI procedure [27]. The momentum 
equations are linked to the continuity equation via 
the method of artificial compressibility. Both 
steady and time-dependent solutions were com- 
puted using an Adams Bashforth-Crank-Nicol- 
son scheme [27]. (For more details see [30].) The 
results presented in this paper were computed 
using 33 X 33 points. This was found to yield 
sufficient resolution for the range of Schmidt 
numbers considered. For the impulse-type residual 
accelerations a variable time step was employed. 
In order to assess the consequences of adopting 
the pseudo-steady state conditions, (9)-(12), we 
varied the aspect ratio of the computational space. 
While, for accelerations between 10“ 4 and 10“ 3 
gal (oriented parallel to the interface) this altered 
the flow pattern in the upper part of the ampoule, 
the computed values of £ differed only by a few 
percent from the values computed with aspect 
ratio 1. For the lower acceleration magnitudes 
(and for orientations perpendicular to the inter- 
face) no significant difference in £ was found. 

PHOENICS embodies a finite volume or finite 
domain formulation [28,31]. It represents the 
governing equations introduced in the previous 
section as a set of algebraic equations. These 
equations represent the consequence of integrating 
the differential equation over the finite volume of 
a computational cell (and, for transient problems, 
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over a finite time) and approximating the resulting 
volume, area and time averages by interpolation. 
For the 2D and axisymmetric calculations dis- 
cussed in this paper we employed a 40 X 39 grid. 
The 3D calculations were performed in a circular 
cylindrical domain with 20 nodes in the radial 
direction, 12 in the azimuthal direction and 39 in 
the axial direction. We found that for the 2D 
calculation less than 40 points in the direction 
parallel to the interface resulted in poor conver- 
gence of the solute field. 

The scheme employed by PHOENICS has the 
same accuracy as a finite difference scheme which 
is of the order Ax, where Ax is the distance 
between the grid nodes. The time scheme is im- 
plicit and thus unconditionally stable. However, 
small time steps are required to obtain accurate 
solutions. For highly non-linear flows the use of 
underrelaxation is necessary to eliminate diver- 
gence and to ensure good convergence of the 
solutions. 

A comparison of the two methods was made 
for different acceleration magnitudes and for dif- 
ferent Schmidt numbers. The spectral method was 
found to cope better with the steep solute gradi- 
ents which occur in the high Schmidt number 
cases. For Sc = 10, and steady acceleration magni- 
tudes of less than 10“ 3 gal, the results of the two 
methods compared favorably. For example, the 
computed concentration values differed by at most 
5%, the difference being less at lower acceleration 
magnitudes. A comparison of interface concentra- 
tion values for accelerations oriented parallel to 
the interface with magnitudes of t/2 X 10“ 6 and 
X 10” 5 , yielded differences of 2% or less. The 
average interface concentration computed by both 
methods being less than 0.1%. The calculated val- 
ues of £ for these two cases were 21.5% and 18.5% 
for the lower acceleration magnitude and 152% 
and 123% for the higher magnitude; in both cases 
the lower value of £ was computed using the 
spectral method. Hence, the percentage difference 
in £ computed by both methods is about a factor 
of 10 (for k = 0.1) higher than the differences in 
the computed values of concentration. 

Steady calculations using the' pseudo-spectral 
method with 33 X 33 collocation points required 
2000-2500 iterations to achieve good convergence 


for the steady calculations. This typically took 
around 5-10 min of (CRAY XMP) CPU time. 
The finite volume calculations were carried out on 
a VAX-785 and required 3-4 h of CPU time. 

4. Results 

4.1. Steady accelerations: 2D 

The temperature field is insensitive to the con- 
vective motion owing to the low Prandtl number 
of the melt (Pr = 10“ 2 ) and the low magnitude of 
the residual accelerations considered. It is char- 
acterized by the conduction profile shown in fig. 
2. Note that there are lateral (radial) as well as 
axial temperature gradients in the system. The 
consequent density gradient is responsible, as we 
shall see, for driving convection in the melt even 
when residual acceleration magnitudes are as low 
as 10“ 6 gal. The temperature fields for all the 
cases listed in table 2 are qualitatively similar, the 
major differences are the magnitudes of the tem- 
perature gradients which reflect the temperature 
difference T H — 7 M . 

The velocity and concentration distributions, 
together with the interfacial £ values, have been 
obtained for various combinations of orientations 


INLET 



Fig. 2. The dimensionless temperature field, 0 , for 2D calcula- 
tions corresponding to operating conditions B-Vl in table 2. 
Note that there is no difference between this and a purely 
conductive profile (i.e. g — 0). 
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Table 3 

Compositional non-uniformity £[%] for Ge:Ga 


Residual 

acceleration 

magnitude 

(gal) 

Orientation 

Ampoule width (cm) 




N 


1 

1 

1 

0.5 

2.0 

Growth rate (p 

ms" 1 ) 




6.5 

3.25 

0.65 

6.5 

6.5 

(A) 1(T 4 

T 

<- 

80 





icr 5 


<— 

92.7 


11.9 

12.0 




X 

70.9 


11.3 





1 

6.4 


0.95 



5X1CT 6 


1 

3.2 







X 

39 







<- 

54.2 





itr 6 



11.3 


2.0 





X 

8.0 





(B) 10“ 4 

T 

1 

36 







t 

32 





]/2 XlO" 5 


X 

110 







<- 

152 (91) 





10" 5 


1 

7.5 

4.6 

0.7 



5^2 X 10" 6 


X 

57 





yT x 10" 6 


X 

10 







t 

4.0 







1 

2.0 







<- 

21.5 (26) 


1.5 



10" 6 


1 

0.7 

0.4 

0.0 


3.8 

10" 7 


T 

1.0 

0.5 

0.2 



(C) 10" 5 

T 

<- 

22.6 




64.5 

10" 6 


«- 

2.3 






e K is the unit vector parallel to g, N is the normal vector to the crystal. 
L=1 cm for all cases. 

Values in parentheses indicate 3D results. 

(A), (B) and (C) refer to the operating conditions listed in table 2. 


and magnitudes of the acceleration vector, growth 
rates, ampoule (02) dimensions and temperature 
differences ?h T m . The results are summarized 
in table 3. 

Figs. 3-5 illustrate the velocity and solute fields 
for three cases where the residual acceleration is 
parallel to the crystal-melt interface and the oper- 
ating conditions correspond to A-Vl in table 2. 
Note that the vectors in these figures are scaled by 
the growth velocity and have been obtained by 
interpolation from the computed values obtained 
on the unevenly spaced Chebyshev mesh. If the 
conditions were purely diffusive, the velocity vec- 
tors would all be the same length and perpendicu- 
lar to the crystal-melt interface. Furthermore, the 


isoconcentrates and crystal-melt interface would 
be parallel. A single roll is calculated at 10 -5 gal, 
resulting in £ = 92.7%. Between 10 -5 and 10 6 gal 
the roll is no longer present and at 10 " 6 gal the 
convective motion results in a barely perceptible 
deflection of the velocity vectors from a purely 
longitudinal flow. The compositional variations in 
the melt are, however, still significant even at 10“ 6 
gal. 

In contrast to the single roll flow depicted in 
fig. 3a, fig. 6 shows that when the same magnitude 
gravity vector is oriented perpendicular to the 
crystal-melt interface there are no rolls. The flow 
is deflected symmetrically about the ampoule axis. 
The associated compositional non-uniformity is 
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g 

INLET 



Fig. 3. (a) The steady flow field produced by a residual acceleration with a magnitude 10 3 gal acting parallel to the crystal-melt 
interface. The maximum speeds are approximately twice the growth speed, (b) The dimensionless solute field, C, associated with the 
flow depicted in (a). The contour interval is 0.9. For this case $ — 92.7%. 

symmetric about the centerline of the ampoule Figs. 7a-7c depict the lateral variation in com- 

and, with £ = 6.4%, is significantly smaller than position of the crystal for gravity orientations 
for the asymmetric case. parallel, at 45 ° and perpendicular to the interface. 



Fig. 4. (a) The steady flow field produced by a residual acceleration with a magnitude 5 X 10~ 6 gal acting parallel to the crystal-melt 
interface. The maximum speeds are slightly greater than the growth speed, (b) The dimensionless solute field, C, associated with the 
flow depicted in (a). The contour interval is 0.7. For this case J = 39%. 







J.I.D. Alexander el al. / Low gravity tolerance of Bridgman— Stockbarger growth. I 


293 



CRYSTAL 



Fig. 5. (a) The steady flow field produced by a residual acceleration with a magnitude 10 -6 gal acting parallel to the crystal-melt 
interface. The maximum speeds are slightly greater than the growth speed, (b) The dimensionless solute field, C, associated with the 

flow depicted in (a). For this case £ = 11.3%. 



CRYSTAL 


Fig. 6. The steady flow field produced by a residual accelera- 
tion with a magnitude ICC 5 gal acting perpendicular to the 
crystal-melt interface, otherwise the operating conditions are 
the same as for figs. 3-5. The contour interval is 0.5. For this 
case | = 6.4%. 


These figures further illustrate the importance of 
ampoule orientation with respect to the residual 
acceleration vector. For a fixed acceleration mag- 
nitude, as the ampoule axis and acceleration vec- 
tor are brought into alignment, the transition from 
asymmetric to symmetric compositional non- 
uniformity is accompanied by a reduction in the 
magnitude of the non-uniformity. For these, and 
other cases corresponding to gallium-doped ger- 
manium crystals, the values of lateral non- 
uniformity £ are listed in table 3 along with the 
associated operating conditions. 

Characterization of melt grown crystals often 
includes radially or laterally averaged crystal com- 
position profiles [32-34], Laterally averaged melt 
composition profiles are shown in fig. 8 for oper- 
ating conditions corresponding to A- VI in table 2 
and for a selection of residual accelerations. Note 
that in fig. 8b these profiles are indistinguishable 
from the pure diffusion case. In addition, these 
profiles, reveal a high concentration gradient re- 
gion adjacent to the interface with an approximate 
width of 0.2 cm. This agrees well with the simple 
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0 0.2 0.4 0.6 0.8 1.0 0 0.2 0.4 0.6 0.8 1.0 0 0.2 0.4 0.6 0.8 1.0 

DISTANCE ALONG INTERFACE (cm) DISTANCE ALONG INTERFACE (cm) DISTANCE ALONG INTERFACE (cm) 


Fig. 7. Lateral variations in melt composition at the interface for operating conditions corresponding to A-Vl in table 2 with 
acceleration magnitudes of 10 -5 , 5X10~ 6 and 10" 6 gal with (a) parallel, (b) 45°, and (c) perpendicular orientation with respect to 

the crystal-melt interface. 

characteristic length obtained from the relation conditions A-V3 (i.e. V M — 0.65 (im s -1 ). A com- 
8 = D/V m . parison of figs. 3a and 9a shows that due to the 

The consequences of reducing the growth rate reduced melt translation rate the buoyancy-driven 
are illustrated by fig. 9 which was obtained for recirculation now dominates. Thus, the axial con- 

1 0 5 gal parallel to the interface and operating centration gradient is reduced (figs. 3b and 9b) 



0 0.2 0.4 0.6 0.8 1.0 0 0.2 0.4 0.6 0.8 1 0 

AXIAL DISTANCE FROM INLET (cm) AXIAL DISTANCE FROM INLET (cm) 

Fig. 8. Laterally averaged, axial composition profiles in the melt for operating conditions corresponding to A-Vl in table 2 for 
accelerations with magnitudes of 10~ 5 , 5X10' 6 and 10 -6 gal oriented (a) parallel and (b) perpendicular to the crystal-melt 
interface. Note that the latter are indistinguishable from the pure diffusion profile. 
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INLET 




Fig. 9. (a) The velocity field associated with a residual acceleration of magnitude 10“ 5 gal acting parallel to the crystal-melt 
interface. The operating conditions correspond to A-V3 in table 2, i.e. the growth rate is 0.65 p m s -1 . (b) The dimensionless solute 
field associated with flow depicted in (a). The contour interval is 0.2. For this case f = 11.9%. 


and a smaller concentration gradient along the 
interface results (in fig. 9b, £ = 11.9% as compared 
to 92.7% in fig. 3b). 

Table 3 also illustrates the effect of changing 
the temperature gradient and the ampoule dimen- 


sions. For an acceleration of lCP 3 gal oriented 
parallel to the interface a reduction of the temper- 
ature difference between the hot zone and the 
crystal-melt interface from 100 to 20 K changes £ 
from 92.7% to 22.6%. An increase in the width of 



Fig. 10. Lateral variations in melt composition at the interface in crystals grown from melts with different Schmidt numbers for 
accelerations with the following magnitudes and orientations with respect to the crystal-melt interface: (a) 10~ 5 gal parallel, (b) 10~ 6 

gal parallel and (c) 10 - 5 gal perpendicular. 
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Fig. 11. Laterally averaged axial composition profiles for melts 
with different Schmidt numbers, for a 10 -4 gal acceleration 
oriented parallel to the crystal-melt interface. 


the ampoule from 1 to 2 cm with the 20 K 
temperature difference yields £ = 64.5%. Clearly, 
for a given steady acceleration the value of £ is 
sensitive to the particular operating conditions 
employed. 

In order to ascertain the importance of the 
diffusivity of the dopant we have also examined 
cases with Sc = 10, 20, 30, and 50 for operating 
conditions otherwise corresponding to A-Vl. (Re- 
call that in the relative weighting of convective 
and diffusive fluxes, an increase in Sc corresponds 
to a reduction in diffusive transport.) The results 
are displayed in fig. 10. It is clear that, all other 
things being equal, the relationship between Sc 
and £ is not linear. In particular, for acceleration 
magnitudes of 10“ 6 and 10“ 5 gal oriented parallel 
to the interface, the maximum value of £ is at- 
tained for Sc = 20, while at 10“ 5 gal oriented 
perpendicular to the interface there is little dif- 
ference in the degree of non-uniformity between 
Sc = 20 and 50. Laterally averaged composition 
profiles for 10“ 4 gal acting parallel to the inter- 
face are depicted for Sc = 10, 20, 30, and 50 in fig. 
11. In comparison to fig. 8 the laterally averaged 
composition profiles at 10“ 4 gal exhibit a flat 


section owing to the increase in convective trans- 
port. 

4.2. Steady accelerations: 3D and axisymmetric 

Five three-dimensional calculations were un- 
dertaken for R = 0.5 cm and V M = 6.5 jiim s“ 1 in 
order to examine the influence of a more realistic 
geometry. In all cases the operating conditions 
corresponded to B-Vl in table 2. Three runs were 
axisymmetric, with g antiparallel and parallel to 
N. These calculations were carried out (for com- 
parison purposes) with || g || = 10“ 4 gal, 10“ 3 gal 
and 10“ 2 gal. Values of the radial segregation 
were found to be within 5% of the values com- 
puted for these conditions by Chang and Brown. 
The compositional non-uniformity £ was found to 
be approximately 10% lower for the 10“ 4 gal 
axisymmetric case for the 2D analog. 

The remaining calculations were fully three-di- 
mensional with g parallel to the crystal-melt in- 
terface. The fully 3D cases were carried out for 
|! g || = X 10“ 5 gal and /2 X 10“ 6 gal. The 
melt isoconcentrates at the interface are shown in 
figs. 12 and 13 for sections cut perpendicular to 
the ampoule axis at z = 1. At the higher value of 
the residual acceleration, £=91%, which is ap- 
proximately half the 2D value. At /2 X 10“ 6 gal, 



Fig. 12. The solute field at the crystal-melt interface conse- 
quent to a </2 xl0~ 5 gal acceleration oriented parallel to the 
interface for the 3D case. The cross section is taken perpendic- 
ular to the ampoule axis. { = 91?c. The contour interval is 0.95. 
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Fig. 13. The solute field at the crystal-melt interface conse- 
quent to yTxlO -6 gal acceleration oriented parallel to the 
interface for the 3D case. The cross section is taken perpendic- 
ular to the ampoule axis. £ <= 26%. The contour interval is 0.26. 

1 = 26%. Thus, for this case, there was little dif- 
ference between the 2D and 3D predictions (see 
table 3). 

4.3. Impulse-type accelerations 

Five cases of impulse-type disturbances were 
examined for a 2D geometry. All impulses were 
superimposed onto a steady background accelera- 
tion which was oriented parallel to the crystal-melt 
interface. The results are summarized in table 4. 


Figs. 14a-14d illustrate the development of the 
solute field following a one second 5 X 10~ 3 gal 
impulse oriented anti-parallel to a background 
acceleration of 10 ~ 6 gal for operating conditions 
corresponding to A-Vl. The effects are long last- 
ing. The velocity field relaxes back to the initial 
state after some 300 s have elapsed. The response 
of the solute field is slower. The effect of the 
impulse is to initially re-orient the flow field. This 
has the effect of first reducing the lateral composi- 
tional nonuniformity. Thirty seconds after the 
termination of the impulse the composition non- 
uniformity is beginning to increase but has the 
opposite sense in comparison to the initial non- 
uniformity (see figs. 14a and 14b). At approxi- 
mately 260 s after the termination of the impulse, 
£ reaches a maximum value of 40%. It than slowly 
decreases in value (figs. 14c and 14d), reverses 
sense and eventually reaches its initial steady level 
(fig. 5b) after more than 2000 s have elapsed. 

A shorter duration (10 _1 s) pulse, with the 
same orientation and magnitude, resulted in a 
maximum deviation of the lateral nonuniformity 
of only 5% from the initial steady level 350 s after 
the termination of the pulse. 

The effects of two one-second pulses separated 
by one second were also calculated. The magni- 
tude of the pulses was 3 X 10 ~ 3 gal. They were 
oriented parallel to the crystal interface and super- 
imposed onto a steady background acceleration of 


Table 4 

Summary of results for impulses 


Magnitude 

(gal) 

Duration 

(s) 

£init 

£ max 

£min 

T v 

(*) 

T ( 

(s) 

T c 

(s) 

5X10“ 3 (ap 10“ 6 ) 

1 

n 

40 

0 

3-400 

260 

2200 

3X10 -3 (ap yT XlO -6 ) 

1 

22 

26 

0 

3-400 

260 

>2000 

3XKT 3 (ap yjl xlO -6 ) 

0.1 

22 

22 

17 

50 

350 

>1000 

3xlO~ 3 (apVTxlO“ 6 ) 

1* 

22 

76 

17 

3-450 

225 

> 2000 

3X10 -3 (ap yT X10 -6 ) 

1 







(VT XKT 6 ) 

No pulse 

22 

24 

22 

300 

100 

>1000 

3 X 10~ 3 (P y/2 X10“ 6 ) 

1 








ap and p denote a pulse anti-parallel and parallel to the background acceleration. All pulses are parallel to the crystal-melt interface. 
The background accelerations are given in parentheses and the asterisk denotes a 1 s pulse that was was repeated once after a 1 s 
interval. 


The relaxation time r v refers to the decay time for the velocity field. t { denotes the time to reach the maximum value of and r c the 
time for the solute field to return to its initial steady state distribution. 
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Fig. 14. The solute field after the termination of a one second duration 0.5X10 2 gal impulse (g ; ) oriented parallel to the 
crystal-melt interface and anti-parallel to a steady 10 -6 gal background: (a) 1 s (£ =11.3%), (b) 30 s (£ = 0.61%), (c) 250 s (£ = 40%) 

and (d) 1250 s (| = 3.4%). 


\Jl X 10" 6 gal oriented parallel to the interface. 
The main effect was to drive £ from 21.5% (the 
initial value) to 76 % after 225 s. 

In addition, double pulses were examined. A 
pulse anti-parallel to the background steady accel- 
eration followed by an equal but opposite pulse 
does not result in a “null” effect. While the flow 
generated by the first pulse is reversed by the 
second pulse there is a net flow following the 
termination of the second pulse. This flow is in the 


same sense as the initial steady flow and results in 
an increase of the lateral non-uniformity in com- 
position to a maximum of 24% at 100 s, whereu- 
pon it decays slowly to its initial value. 

5. Summary and discussion 

The effects on dopant uniformity during direc- 
tional solidification of steady and impulse-type 
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residual accelerations have been calculated using 
2D and 3D “generic” models of the Bridgman- 
Stockbarger crystal growth technique. The salient 
results of our calculations for materials with prop- 
erties and growth conditions corresponding to 
those listed in table 2 can be summarized as 
follows: 

For a fixed growth rate, the amount of lateral 
non-uniformity in composition is very sensitive to 
the orientation of the steady component of the 
residual gravity vector. The worst case appears to 
be when the acceleration vector is parallel to the 
crystal interface. At growth rates on the order of 
microns per second, this orientation can lead to 
non-uniformities of 10-20% when the magnitude 
of the acceleration is 10“ 6 gal. If, however the 
growth rate is lowered by an order of magnitude, 
the non-uniformity is reduced significantly. A 
steady background level on the order of 10“ 6 gal 
can be tolerated for a wide range of the operating 
conditions examined provided that the accelera- 
tion vector is aligned with the axis of the growth 
ampoule. 

When g is oriented parallel to the crystal-melt 
interface, the maximum compositional non-uni- 
formity occurs between 10“ 4 and 10“ 5 gal com- 
pared to 10“ 2 gal for the axisymmetric case. It is 
evident upon examination of table 3 and figs. 3-7 
that for a given orientation the relationship be- 
tween the magnitude of the gravity vector and the 
amount of compositional non-uniformity is clearly 
not linear. (This result was first obtained by Chang 
and Brown [13] for the axisymmetric case.) 

The effect of a reduction in growth rate on the 
velocity profile is to shift the center of the convec- 
tive roll (see figs. 9a and 3a). This can be explained 
in terms of the superposition of the uniform solid- 
ification flow and a single convective roll. When 
convective velocities have the same order of mag- 
nitude as the solidification velocity the roll center 
is shifted to one side because only here does the 
addition of the convective and uniform solidifica- 
tion flow result in no motion. 

For the cases examined, a reduction in com- 
positional non-uniformity was also associated with 
a reduction in t h t m . The presence of the 
adiabatic zone results in lateral as well as axial 
temperature gradients. Thus, convection always 


occurs in this system since there is always a loca- 
tion where the temperature gradient has a compo- 
nent normal to the acceleration vector. We com- 
pared our results with those of Chang and Brown’s 
axisymmetric calculations [13]. Our 3D axisym- 
metric computations were found to be in agree- 
ment with their work. In addition we carried out 
full 3D ( non-axisymmetric ) steady calculations in 
order to calibrate our 2D results. For one set of 
examples we found that with fl X 10“ 5 gal paral- 
lel to the interface the percentage compositional 
non-uniformity predicted by the 2D calculation 
(I = 152%) was 50% higher than predicted by the 
full 3D calculation. This difference may be attri- 
buted to the increased surface to volume ratio 
which increases the effect of the “no-slip” 
boundary condition. The presence of the rigid 
walls somewhat retards the flow and, thus, in 
comparison to the 2D case, the degree of solute 
redistribution is decreased. A comparison between 
2D and 3D calculations for -Jl X 10“ 6 gal re- 
vealed no significant difference in £. At this mag- 
nitude of the gravity vector, convection is local- 
ized and weak. The increase in surface to volume 
ratio for this 3D case thus has little influence on 
the transport conditions. 

There are some instances in the literature where 
radially or laterally averaged crystal composition 
profiles are used to compare space and earth 
grown samples [32-33]. In their interpretation of 
laterally averaged composition profiles in earth- 
and space-grown Ge : Ga crystals [32], Witt et al. 
concluded that for the space-grown crystal, con- 
vection was absent during growth. The theory of 
Tiller et al. [34], which considers only diffusion, 
gives a characteristic distance z* = D/kV at 
which the crystal composition reaches a value 
given by c* = c 0 [l — (1 — k)e~ 1 ]. The distance at 
which this composition was reached in the earth- 
and space-grown samples was compared with the 
theoretical distance z*. Given the uncertainty as- 
sociated with the relevant physical properties and 
the fact that the growth rate was not constant (as 
assumed in the theory of Tiller et al.), a range was 
given for the theoretical distance z * . It was found 
that for the space-grown sample the measured 
characteristic distance fell within this range. The 
composition c * was not reached in the earth-grown 
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sample. Thus, it seems appropriate to argue that 
convective transport is absent if the measured 
characteristic distance is equal to z*. However, an 
analogous argument applied to the averaged melt 
composition profiles of figs. 8 and 11 leads us to a 
different conclusion. It has been proposed [35,36] 
that an appropriate measure of the influence of 
convective transport is the ratio, A, of the solute 
boundary layer thickness 8 to the characteristic 
distance D/V u . Corresponding transport regimes 
are then defined as [35,36]: 

4«1: convective, 

A = 1 : convective-diffusive, 

A > 1 : diffusive. 

As a working definition we take 5 to be the 
distance from the interface where the rate of 
change of concentration gradient becomes small 
and evaluate A for the profiles given in figs. 8 and 
11. In fig. 11, the approximate values of A range 
from 1 (Sc = 10) to 0.75 (Sc = 40). For all the 
profiles in figs. 8, which are virtually indis- 
tinguishable from the purely diffusive profiles, A 
is greater than 1. However, the profiles shown in 
fig. 8a are associated with significant lateral com- 
positional non-uniformity that arises due to con- 
vection. Thus, even when the laterally averaged 
melt composition profiles have so-called diffusive 
(longitudinal) profiles, it may not be safe to as- 
sume that convection was absent. Thus, we pro- 
pose that radially or laterally averaged composi- 
tion profiles are alone insufficient to describe the 
extent of residual convection during solidification 
in a spacecraft environment. 

The effects of impulse-type disturbances can be 
severe and can extend for a long time (on the 
order of 10 3 s) after the termination of the im- 
pulse. This is not surprising if one considers char- 
acteristic times for momentum and solute diffu- 
sion given by t m = L 2 /v and t s = L 2 /D, respec- 
tively. For L — 0.5 cm (corresponding to the maxi- 
mum distance from a wall in the 2D system con- 
sidered) these characteristic times are t m ~ 200 s, 
and r s ~ 2000 s, which compare well with the 
order of magnitude of the decay times given in 
table 4. At growth rates of 6.5 pm s -1 , a pulse 
with a second duration, or a combination of such 
pulses has a drastic effect on the segregation levels 


at pulse amplitudes of 10 ~ 3 gal. The nature of the 
response depends on the magnitude, direction and 
duration of the impulse, and whether sequential, 
opposing impulses are involved. A so-called “com- 
pensating” double pulse will not result in com- 
pletely offsetting effects. For the case we ex- 
amined, however, the resulting compositional 
non-uniformity was not as severe as for sequential 
pulses with the same orientation. It should also be 
considered that the response will depend on the 
particular operating conditions employed. For ex- 
ample, the response to steady acceleration was 
shown to vary with the growth rate, and a similar 
dependence may exist for impulses. This remains 
to be verified. Further investigation of more realis- 
tic impulses (g-jitter) is necessary since the re- 
sponse of the system appears to depend on the 
nature of the impulse and our results indicate that 
impulses appear to have important consequences 
for transient behavior in crystal growth systems. 

In addition to the results obtained for Ge : Ga, 
we have also examined the effects of different 
physical properties by examining a range of 
Schmidt numbers. Our results show that, for a 
given steady acceleration the sensitivity of the 
process also depends on the operating conditions 
and the physical properties of the melt. That is, 
the sensitivity to the residual acceleration depends 
on temperature conditions, the growth rate, and 
the ratio of the kinematic viscosity to the solute 
diffusivity. The fact that in a spacecraft buoyancy 
forces may be reduced by 5-6 orders of magni- 
tude does not imply that the environment may be 
considered purely diffusive. Careful evaluation of 
the interaction between the operating conditions, 
physical properties and spacecraft acceleration en- 
vironment must be made in order to ensure opti- 
mum use of limited spaceflight opportunities. 

Space flight experiments involving germanium 
have been carried out on Skylab [32,33,37] and the 
Apollo-Soyuz Test Project (ASTP) [38,39], The 
latter involved the directional solidification of 
seeded Ge melts. Melts were doped with galium 
and some were doped with 1% Si and 0.001% Sb 
[39]. These results have generated some con- 
troversy [1]; lateral compositional profiles re- 
vealed strong asymmetric non-uniformities in the 
space grown crystals. Lateral variations were also 
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observed in samples grown under terrestrial condi- 
tions but were much less pronounced [1]. In order 
to explain these observations it has been argued 
that the asymmetric redistribution of the dopant 
can be ascribed to “barometric diffusion” of the 
solutes due to the acceleration gradient in the melt 
arising from the rotational motion of the spacecraft 
[38,39]. The basis of the argument, however, ap- 
pears to ignore the presence of gravity gradient 
and atmospheric drag effects and does not ex- 
plicitly account for the spacecraft attitude mo- 
tions. If it is assumed that during this phase of the 
ASTP mission the spacecraft flew in a radial atti- 
tude [9] (i.e. it continuously rotated relative to a 
fixed geocentric frame), then a steady residual 
acceleration vector would result from the gravity 
gradient, atmospheric drag and centrifugal contri- 
butions to the relative acceleration. It has been 
demonstrated in this work that melt convection 
can occur in response to a steady residual acceler- 
ation of rather small magnitudes. Furthermore, 
whenever the acceleration vector is not aligned 
with the ampoule axis strong asymmetries in the 
solute composition profiles can occur. Since the 
dependence of composition uniformity on convec- 
tion is not linear, and since small asymmetries in 
temperature profiles can lead to asymmetric con- 
vection [40], there is always the possibility that the 
observed non-uniformities in the terrestrially 
grown samples will be smaller than those observed 
in space grown samples. Thus, the probability that 
the strong asymmetries in composition observed 
in the space-grown samples can be attributed to 
residual acceleration effects should not be dis- 
counted. 

It should be emphasized that while some of the 
conditions examined have shown great sensitivity 
to residual acceleration, we have also shown that 
there are combinations of operating conditions for 
which the effects of realistic residual acceleration 
levels are small. It is therefore clear that residual 
gravity effects should be seriously considered when 
planning optimal operating conditions for space 
experiments. Furthermore, since our results are 
generic in nature they cannot be substituted for a 
detailed numerical examination of a particular 
Bridgman-type crystal growth system, although 
they can be expected to indicate trends in the 


expected response to the steady residual accelera- 
tion environment. As far as the response to time- 
dependent accelerations is concerned, impulses 
have been shown to be important, but more work 
is needed in this area, especially in terms of a 
more realistic representation of the residual accel- 
eration vector. Continuous periodic accelerations 
are also a concern [21], particularly in the low 
frequency range (~ 10 ~ 3 Hz) and will be examined 
in a companion paper [41]. 
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